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

Скрипт Open 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

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)}$ (Раздел See Van Loan 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$ комплексный вектор, вертят факторы. Вертеть факторы, $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.

Как замечено в графике ниже, ошибка является в допуске 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

Осведомленная о типе функция битного реверсирования

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, чтобы Идентифицировать Переполнение

Чтобы оснастить код 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 / макс. значения.

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], [валлийские 1969]. Однако, если вы масштабируете вход к каждому из этапов БПФ 1/2, можно получить полное масштабирование 1/n и произвести отношение сигнал-шум, пропорциональное n [Oppenheim & Schafer 1989, уравнение 9.105], [валлийские 1969].

Эффективным способом масштабироваться 1/2 в фиксированной точке являются к сдвигу вправо данные. Для этого вы используете функцию арифметики права сдвига разряда 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'});

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

Очистка

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

fipref(FIPREF_STATE);

Ссылки

Ссуда Чарльза Вана, вычислительные среды для быстрого преобразования Фурье, SIAM, 1992.

Клив Moler, числовое вычисление с MATLAB, SIAM, 2004, анализ Фурье главы 8.

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

Питер Д. Уэлч, "Анализ ошибок Быстрого преобразования Фурье Фиксированной точки", IEEE® Transactions на Аудио и Электроакустике, издании AU-17, № 2, июнь 1969, стр 151-157.

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