Алгоритм FFT учебника

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

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

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

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

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

В этом примере используется метод прореживания за единицу времени, показанный в алгоритме 1.6.2 на стр. 45 книги «Вычислительные структуры для быстрого преобразования Фурье» Чарльза Ван Займа.

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

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

$\mathit{w}{\mathit{=}}_{\mathit{}}^{\mathrm{wn\; (}}$длинный) (См. раздел 1.4.11 «Кредит фургона»)

для $\mathit{q}\mathrm{=}\mathit{1}$: t

для $\mathit{k}\mathrm{=}\mathit{0}:\mathrm{}$r-1

для $\mathit{j}\mathrm{=}{\mathit{0}}_{:}\mathrm{L}$* -1

конец

конец

конец

Алгоритм учебника использует индексацию на основе нуля. ${\mathit{}}_{\mathit{Fn}}$ является матрицей n-на-n преобразования Фурье, ${\mathit{}}_{\mathit{Pn}}$ является матрицей перестановки n-на-n битов-реверсирования, и $\mathit{w}$ является комплексным вектором множителей твиддла. Множители widdle, $\mathit{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')

Проверка кода с плавающей запятой

Для реализации алгоритма в 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.

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

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'});

Преобразование функций для использования таблиц типов

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

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

Таблица исходных типов

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

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

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

Функция реверсирования битов с учетом типа

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

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'});

Создание таблицы типов фиксированных точек

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

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 Instrumentation для определения переполнений

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

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

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

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 для записи значений min/max.

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 в фиксированной точке является смещение данных вправо. Для этого используется арифметическая функция «bit shift right» bitsra. После масштабирования каждого этапа БПФ и оптимизации вычисления индексной переменной алгоритм становится следующим:

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'});

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

Очистить

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

fipref(FIPREF_STATE);

Ссылки

Charles Van Loan, Вычислительные рамки для быстрого преобразования Фурье, SIAM, 1992.

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

Алан В. Оппенгейм и Рональд В. Шефер, дискретная обработка временных сигналов, Прентис Холл, 1989.

Питер Д. Уэлч, «Анализ ошибок быстрого преобразования Фурье с фиксированной точкой», IEEE ® Transactions on Audio and Electroacustics, том AU-17, № 2, июнь 1969, стр. 151-157.