Преобразуйте быстрое преобразование Фурье (FFT) в фиксированную точку

Открыть Live Script

В этом примере показано, как преобразовать версию учебника алгоритма быстрого преобразования Фурье (FFT) в код MATLAB ® с фиксированной точкой.

Запустите следующий код, чтобы захватить текущие состояния и сбросить глобальные состояния.

FIPREF_STATE = get(fipref);
reset(fipref)

БПФ учебника

БПФ является комплексным линейным преобразованием из временного интервала в частотный диапазон. Для примера, если создать вектор как сумму двух синусоид и преобразовать его с помощью БПФ, можно увидеть peaks частот на графике БПФ величины.

n = 64;                                     % Number of points
Fs = 4;                                     % Sampling frequency in Hz
t  = (0:(n-1))/Fs;                          % Time vector
f  = linspace(0,Fs,n);                      % Frequency vector
f0 = 0.2; f1 = 0.5;                         % Frequencies, in Hz
x0 = cos(2*pi*f0*t) + 0.55*cos(2*pi*f1*t);  % Time-domain signal
x0 = complex(x0);                           % The textbook algorithm requires the input to be complex
y0  = fft(x0);                              % Frequency-domain transformation fft() is a MATLAB built-in function

fi_fft_demo_ini_plot(t,x0,f,y0);            % Plot the results from fft and time-domain signal

Figure contains 2 axes. Axes 1 contains an object of type line. This object represents x0. Axes 2 contains an object of type line. This object represents abs(fft(x0)).

Peaks на частотном графике 0,2 и 0,5 Гц соответствуют двум синусоидам сигнала временной области на этих частотах.

Обратите внимание на отражённый peaks с частотой 3,5 и 3,8 Гц. Когда вход к БПФ является реальным, как это в этом случае, то выход $y$ является сопряженным-симметричным:

$y (k) = conj (y (n - k))$

Существует много различных реализаций БПФ, каждая из которых имеет свои затраты и преимущества. Вы можете обнаружить, что другой алгоритм лучше для вашего приложения, чем тот, который приведен здесь. Этот алгоритм предоставляет вам пример того, как можно начать собственное исследование.

Этот пример использует децимацию во времени с единичным шагом БПФ показанную в Алгоритме 1.6.2 на странице 45 книги «Вычислительные Среды для быстрого преобразования Фурье» Чарльза Ван Лоуна.

В псевдокоде алгоритм в учебнике следующий:

Алгоритм 1.6.2. Если $x$ является комплексным вектором длины $n$ и $n = 2^{t}$ , затем следующий алгоритм перезаписывает $x$ с $F_{n} x$ .

$x = P_{n} x$

$w = w_{n}^{(long)}$ (См. Раздел 1.4.11 «Заем фургона».)

для $q = 1 : t$

$L = 2^{q}; r = n / L; L_{*} = L / 2;$

для $k = 0 : r - 1$

для $j = 0 : L_{*} - 1$

$τ = w (L_{*} - 1 + j) \cdot x (kL + j + L_{*})$

$x (kL + j + L_{*}) = x (kL + j) - τ$

$x (kL + j) = x (kL + j) + τ$

конец

Алгоритм учебника использует нулевую индексацию. $F_{n}$ является n-на-n матрицей Преобразования Фурье, $P_{n}$ является n-на-n матрицей сочетания обратного бита, и $w$ является комплексным вектором twiddle множителей. Коэффициенты твиддла, $w$ , являются комплексными корнями единицы, вычисленными следующим алгоритмом:

function w = fi_radix2twiddles(n) 
%FI_RADIX2TWIDDLES  Twiddle factors for radix-2 FFT example.
%   W = FI_RADIX2TWIDDLES(N) computes the length N-1 vector W of
%   twiddle factors to be used in the FI_M_RADIX2FFT example code.
%
%   See also FI_RADIX2FFT_DEMO.

%   Reference:
%
%   Twiddle factors for Algorithm 1.6.2, p. 45, Charles Van Loan,
%   Computational Frameworks for the Fast Fourier Transform, SIAM,
%   Philadelphia, 1992.
%
%   Copyright 2003-2011 The MathWorks, Inc.
%     

t = log2(n);
if floor(t) ~= t
  error('N must be an exact power of two.');
end

w = zeros(n-1,1);
k = 1;
L = 2;
% Equation 1.4.11, p. 34
while L <= n
  theta = 2*pi/L;
  % Algorithm 1.4.1, p. 23
  for j = 0:(L/2 - 1)
    w(k) = complex(cos(j*theta),-sin(j*theta));
    k = k + 1;
  end
  L = L*2;
end

figure(gcf)
clf
w0 = fi_radix2twiddles(n);
polarplot(angle(w0),abs(w0),'o')
title('Twiddle Factors: Complex roots of unity')

Figure contains an axes. The axes with title Twiddle Factors: Complex roots of unity contains an object of type line.

Проверьте код с плавающей точкой

Для реализации алгоритма в MATLAB ® можно использовать fi_bitreverse функция для обратного бита последовательности входа. Необходимо добавить единицу к индексам, чтобы преобразовать их из основанных на нуле в основанные на одном.

function x = fi_m_radix2fft_algorithm1_6_2(x, w)
%FI_M_RADIX2FFT_ALGORITHM1_6_2  Radix-2 FFT example.
%   Y = FI_M_RADIX2FFT_ALGORITHM1_6_2(X, W) computes the radix-2 FFT of
%   input vector X with twiddle-factors W.  Input X is assumed to be
%   complex.
%
%   The length of vector X must be an exact power of two.
%   Twiddle-factors W are computed via
%      W = fi_radix2twiddles(N)
%   where N = length(X).
%
%   This version of the algorithm has no scaling before the stages.
%
%   See also FI_RADIX2FFT_DEMO, FI_M_RADIX2FFT_WITHSCALING.

%   Reference:
%     Charles Van Loan, Computational Frameworks for the Fast Fourier
%     Transform, SIAM, Philadelphia, 1992, Algorithm 1.6.2, p. 45.
%
%   Copyright 2004-2015 The MathWorks, Inc.

    n = length(x);  t = log2(n);
    x = fi_bitreverse(x,n);
    for q=1:t
        L = 2^q; r = n/L; L2 = L/2;
        for k = 0:(r-1)
            for j = 0:(L2-1)
                temp          = w(L2-1+j+1) * x(k*L+j+L2+1);
                x(k*L+j+L2+1) = x(k*L+j+1)  - temp;
                x(k*L+j+1)    = x(k*L+j+1)  + temp;
            end
        end
    end
end

Визуализация

Чтобы убедиться, что вы правильно реализовали алгоритм в MATLAB ®, запустите через него известный сигнал и сравните результаты с результатами, полученными функцией MATLAB ® FFT.

Как видно на графике ниже, ошибка находится в пределах допуска встроенной функции MATLAB ® FFT, проверяя, что вы правильно реализовали алгоритм.

y = fi_m_radix2fft_algorithm1_6_2(x0, w0);

fi_fft_demo_plot(real(x0),y,y0,Fs,'Double data', ...
    {'FFT Algorithm 1.6.2','Built-in FFT'});

Преобразуйте функции в таблицы типов

Чтобы отделить типы данных от алгоритма:

Составьте таблицу определений типов данных.
Измените код алгоритма, чтобы использовать типы данных из этой таблицы.

Этот пример показывает итерационные шаги путем создания различных файлов. На практике можно внести итерационные изменения в тот же файл.

Таблица типов

Составьте таблицу типов с помощью структуры с прототипами для переменных, установленными на их исходные типы. Используйте типы опорной структуры, чтобы подтвердить, что начальное преобразование выполнено правильно, и программно переключить функцию между типами с плавающей и фиксированной точками. Переменные индекса автоматически преобразуются в целые числа с помощью Coder™ MATLAB ®, поэтому вам не нужно указывать их типы в таблице.

Задайте значения прототипа как пустые ([]), так как используются типы данных, но не значения.

function T = fi_m_radix2fft_original_types()
%FI_M_RADIX2FFT_ORIGINAL_TYPES Types Table Example

%   Copyright 2015 The MathWorks, Inc.

    T.x = double([]);
    T.w = double([]);
    T.n = double([]);

end

Функция алгоритма с поддержкой типа

Добавьте T таблицы типов как вход в функцию и использовать ее, чтобы привести переменные к конкретному типу, сохраняя при этом тело алгоритма неизменным.

function x = fi_m_radix2fft_algorithm1_6_2_typed(x, w, T)
%FI_M_RADIX2FFT_ORIGINAL_TYPED Radix-2 FFT example.
%   Y = FI_M_RADIX2FFT_ALGORITHM1_6_2_TYPED(X, W, T) computes the radix-2
%   FFT of input vector X with twiddle-factors W.  Input X is assumed to be
%   complex.
%
%   The length of vector X must be an exact power of two.
%   Twiddle-factors W are computed via
%      W = fi_radix2twiddles(N)
%   where N = length(X).
%
%   T is a types table to cast variables to a particular type, while keeping
%   the body of the algorithm unchanged.
%
%   This version of the algorithm has no scaling before the stages.
%
%   See also FI_RADIX2FFT_DEMO, FI_M_RADIX2FFT_WITHSCALING.
%
%   Reference:
%     Charles Van Loan, Computational Frameworks for the Fast Fourier
%     Transform, SIAM, Philadelphia, 1992, Algorithm 1.6.2, p. 45.

%   Copyright 2015 The MathWorks, Inc.
%
%#codegen

    n = length(x);
    t = log2(n);
    x = fi_bitreverse_typed(x,n,T);
    LL = cast(2.^(1:t),'like',T.n);
    rr = cast(n./LL,'like',T.n);
    LL2 = cast(LL./2,'like',T.n);
    for q=1:t
        L = LL(q);
        r = rr(q);
        L2 = LL2(q);
        for k = 0:(r-1)
            for j = 0:(L2-1)
                temp          = w(L2-1+j+1) * x(k*L+j+L2+1);
                x(k*L+j+L2+1) = x(k*L+j+1)  - temp;
                x(k*L+j+1)    = x(k*L+j+1)  + temp;
            end
        end
    end
end

Функция bit-reversal с поддержкой типа

function x = fi_bitreverse_typed(x,n0,T)
%FI_BITREVERSE_TYPED  Bit-reverse the input.
%   X = FI_BITREVERSE_TYPED(x,n,T) bit-reverse the input sequence X, where
%   N=length(X).
%
%   T is a types table to cast variables to a particular type, while keeping
%   the body of the algorithm unchanged.
%
%   See also FI_RADIX2FFT_DEMO.

%   Copyright 2004-2015 The MathWorks, Inc.
%
%#codegen
n = cast(n0,'like',T.n);
nv2 = bitsra(n,1);
j = cast(1,'like',T.n);
for i = 1:(n-1)
  if i < j
    temp = x(j);
    x(j) = x(i);
    x(i) = temp;
  end
  k = nv2;
  while k < j
    j(:) = j-k;
    k = bitsra(k,1);
  end
  j(:) = j+k;
end

Проверьте измененную функцию

Каждый раз, когда вы изменяете свою функцию, проверяйте, что результаты все еще соответствуют вашему базовому уровню. Поскольку исходные типы использовались в таблице типов, выходы должны быть идентичными. Это подтверждает, что вы произвели преобразование, чтобы правильно отделить типы от алгоритма.

T1 = fi_m_radix2fft_original_types(); % Get original data types declared in table

x = cast(x0,'like',T1.x);
w = cast(w0,'like',T1.w);

y = fi_m_radix2fft_algorithm1_6_2_typed(x, w, T1);

fi_fft_demo_plot(real(x),y,y0,Fs,'Double data', ...
    {'FFT Algorithm 1.6.2','Built-in FFT'});

Figure contains 3 axes. Axes 1 contains an object of type line. This object represents Double data. Axes 2 contains 2 objects of type line. These objects represent FFT Algorithm 1.6.2, Built-in FFT. Axes 3 contains an object of type line. This object represents abs(error).

Составьте таблицу Фиксированных точек

Составьте таблицу фиксированных точек с помощью структуры с прототипами для переменных. Задайте значения прототипа как пустые ([]), так как используются типы данных, но не значения.

function T = fi_m_radix2fft_fixed_types()
%FI_M_RADIX2FFT_FIXED_TYPES Example function

%   Copyright 2015 The MathWorks, Inc.

    T.x = fi([],1,16,14); % Picked the following types to ensure that the
    T.w = fi([],1,16,14); % inputs have maximum precision and will not
                          % overflow
    T.n = int32([]);      % Picked int32 as n is an index

end

Идентифицируйте проблемы с фиксированной точкой

Теперь попробуйте преобразовать входные данные в фиксированную точку и посмотрите, выглядит ли алгоритм все еще хорошо. В этом первом проходе вы используете все значения по умолчанию для подписанных данных с фиксированной точкой при помощи fi конструктор.

T2 = fi_m_radix2fft_fixed_types(); % Get fixed point data types declared in table

x = cast(x0,'like',T2.x);
w = cast(w0,'like',T2.w);

Повторно запустите тот же алгоритм с входами с фиксированной точкой.

y  = fi_m_radix2fft_algorithm1_6_2_typed(x,w,T2);
fi_fft_demo_plot(real(x),y,y0,Fs,'Fixed-point data', ...
    {'Fixed-point FFT Algorithm 1.6.2','Built-in FFT'});

Обратите внимание, что график величины (центр) БПФ с фиксированной точкой не напоминает график встроенного БПФ. Ошибка (нижний график) намного больше, чем вы ожидаете увидеть для ошибки округления, поэтому вероятно, что произошло переполнение.

Используйте Min/Max Инструментирования, чтобы идентифицировать переполнения

Чтобы инструментализировать код MATLAB ®, создайте MEX-функцию из функции MATLAB ® с помощью buildInstrumentedMex команда. Входы для buildInstrumentedMex являются такими же, как и входы в fiaccel, но buildInstrumentedMex не имеет fi ограничения по объектам. Область выхода buildInstrumentedMex является MEX-функцией с вставленным инструментированием, поэтому при запуске MEX-функции моделируемые минимальное и максимальное значения записываются для всех именованных переменных и промежуточных значений.

The '-o' опция используется для присвоения имени MEX-функции. Если на '-o' опция не используется, тогда MEX-функция является именем функции MATLAB ® с '_mex' приложенный. Можно также назвать MEX-функцию так же, как и функцию MATLAB ®, но необходимо помнить, что MEX-функции имеют приоритет над функциями MATLAB ®, и поэтому изменения функции MATLAB ® не будут запускаться до тех пор, пока MEX-функция не будет регенерирована, или MEX-функция не будет удалена и очищена.

Создайте вход с масштабированным типом данных double, чтобы его значения достигли полной области значений, и вы можете идентифицировать потенциальные переполнения.

function T = fi_m_radix2fft_scaled_fixed_types()
%FI_M_RADIX2FFT_SCALED_FIXED_TYPES Example function

%   Copyright 2015 The MathWorks, Inc.

    DT = 'ScaledDouble';                % Data type to be used for fi
                                        % constructor
    T.x = fi([],1,16,14,'DataType',DT); % Picked the following types to
    T.w = fi([],1,16,14,'DataType',DT); % ensure that the inputs have
                                        % maximum precision and will not
                                        % overflow
    T.n = int32([]);                    % Picked int32 as n is an index

end

T3 = fi_m_radix2fft_scaled_fixed_types(); % Get fixed point data types declared in table

x_scaled_double = cast(x0,'like',T3.x);
w_scaled_double = cast(w0,'like',T3.w);

buildInstrumentedMex fi_m_radix2fft_algorithm1_6_2_typed ...
    -o fft_instrumented -args {x_scaled_double w_scaled_double T3}

Запустите инструментальную MEX-функцию, чтобы записать минимальные/максимальные значения.

y_scaled_double = fft_instrumented(x_scaled_double,w_scaled_double,T3);

Отображение результатов инструментирования.

showInstrumentationResults fft_instrumented

Из результатов инструментирования видно, что при присвоении переменной были переполнения x.

Измените алгоритм, чтобы решить проблемы с фиксированной точкой

Величина отдельного интервала в БПФ увеличивается, самое большее, в n раз, где n - длина БПФ. Следовательно, масштабируя ваши данные на 1/n, можно предотвратить переполнение для любого входа. Когда вы масштабируете только вход в первый каскад БПФ длины-n на 1/n, вы получаете отношение сигнал/шум, пропорциональное n ^ 2 [Oppenheim & Schafer 1989, уравнение 9.101], [Welch 1969]. Однако, если вы масштабируете вход в каждый из каскадов БПФ на 1/2, можно получить общее масштабирование 1/n и получить отношение сигнал/шум, пропорциональное n [Oppenheim & Schafer 1989, уравнение 9.105], [Welch 1969].

Эффективный способ масштабирования на 1/2 в фиксированной точке - сдвинуть данные вправо. Для этого вы используете арифметическую функцию сдвиг right bitsra. После масштабирования каждого этапа FFT и оптимизации расчета индексной переменной, ваш алгоритм становится:

function x = fi_m_radix2fft_withscaling_typed(x, w, T)
%FI_M_RADIX2FFT_WITHSCALING  Radix-2 FFT example with scaling at each stage.
%   Y = FI_M_RADIX2FFT_WITHSCALING_TYPED(X, W, T) computes the radix-2 FFT of
%   input vector X with twiddle-factors W with scaling by 1/2 at each stage.
%   Input X is assumed to be complex.
%
%   The length of vector X must be an exact power of two.
%   Twiddle-factors W are computed via
%      W = fi_radix2twiddles(N)
%   where N = length(X).
%
%   T is a types table to cast variables to a particular type, while keeping
%   the body of the algorithm unchanged.
%
%   This version of the algorithm has no scaling before the stages.
%
%   See also FI_RADIX2FFT_DEMO.

%   Reference:
%     Charles Van Loan, Computational Frameworks for the Fast Fourier
%     Transform, SIAM, Philadelphia, 1992, Algorithm 1.6.2, p. 45.
%
%   Copyright 2004-2015 The MathWorks, Inc.
%
%#codegen

    n = length(x);  t = log2(n);
    x = fi_bitreverse(x,n);

    % Generate index variables as integer constants so they are not computed in
    % the loop.
    LL = cast(2.^(1:t),'like',T.n);
    rr = cast(n./LL,'like',T.n);
    LL2 = cast(LL./2,'like',T.n);
    for q = 1:t
        L = LL(q); r = rr(q); L2 = LL2(q);
        for k = 0:(r-1)
            for j = 0:(L2-1)
                temp          = w(L2-1+j+1) * x(k*L+j+L2+1);
                x(k*L+j+L2+1) = bitsra(x(k*L+j+1) - temp, 1);
                x(k*L+j+1)    = bitsra(x(k*L+j+1) + temp, 1);
            end
        end
    end
end

Запустите масштабированный алгоритм с данными с фиксированной точкой.

x = cast(x0,'like',T3.x);
w = cast(w0,'like',T3.w);

y = fi_m_radix2fft_withscaling_typed(x,w,T3);

fi_fft_demo_plot(real(x), y, y0/n, Fs, 'Fixed-point data', ...
    {'Fixed-point FFT with scaling','Built-in FFT'});

Figure contains 3 axes. Axes 1 contains an object of type line. This object represents Fixed-point data. Axes 2 contains 2 objects of type line. These objects represent Fixed-point FFT Algorithm 1.6.2, Built-in FFT. Axes 3 contains an object of type line. This object represents abs(error).

Figure contains 3 axes. Axes 1 contains an object of type line. This object represents Fixed-point data. Axes 2 contains 2 objects of type line. These objects represent Fixed-point FFT with scaling, Built-in FFT. Axes 3 contains an object of type line. This object represents abs(error).

Вы можете увидеть, что масштабированный алгоритм БПФ с фиксированной точкой теперь соответствует встроенному БПФ с допуском, который ожидается для 16-битных данных с фиксированной точкой.

Очистка

Запустите следующий код, чтобы восстановить глобальные состояния.

fipref(FIPREF_STATE);

Ссылки

Charles Van Loan, Computational Framework for the Быстрое Преобразование Фурье, SIAM, 1992.

Cleve Moler, Численные вычисления с MATLAB, SIAM, 2004, Глава 8 Анализ Фурье.

Alan V. Oppenheim and Ronald W. Schafer, Discrete-Time Signal Processing, Prentice Hall, 1989.

Peter D. Welch, «A Fixed-Point Fast Fourier Transform Error Analysis», IEEE ® Транзакции по аудио и электроакустике, том AU-17, № 2, июнь 1969, стр. 151-157.

Документация