FFT алгоритм в MATLAB

Question

Я ищу реализовать алгоритм БПФ на микроконтроллерах, поэтому я хочу, чтобы имитировать коды, прежде чем на самом деле использовать его

Я получил 2 примера, которые я преобразованы в MatLab коды, но результат просто не то, что я ожидал

Вот коды:

function [ H ] = fft_2(g) 
%FFT2 Summary of this function goes here 
% Detailed explanation goes here 

NUMDATA = length(g); 
NUMPOINTS = NUMDATA/2; 
N = NUMPOINTS; 
% for(k=0; k<N; k++) 
% { 
% IA[k].imag = -(short)(16383.0*(-cos(2*pi/(double)(2*N)*(double)k))); 
% IA[k].real = (short)(16383.0*(1.0 - sin(2*pi/(double)(2*N)*(double)k))); 
% IB[k].imag = -(short)(16383.0*(cos(2*pi/(double)(2*N)*(double)k))); 
% IB[k].real = (short)(16383.0*(1.0 + sin(2*pi/(double)(2*N)*(double)k))); 
% } 
for k=0:(N-1) 
    IA(k+1,2) = -floor(16383.0*(-cos(2*pi/(2*N)*k))); 
    IA(k+1,1) = floor(16383.0*(1.0 - sin(2*pi/(2*N)*k))); 
    IB(k+1,2) = -floor(16383.0*(cos(2*pi/(2*N)*k))); 
    IB(k+1,1) = floor(16383.0*(1.0 + sin(2*pi/(2*N)*k))); 
end 
% Note, IA(k) is the complex conjugate of A(k) and IB(k) is the complex conjugate of 
% B(k). 
% *********************************************************************************/ 
% #include <math.h> 
% #include ”params1.h” 
% #include ”params.h” 
% extern short g[]; 
% void dft(int, COMPLEX *); 
% void split(int, COMPLEX *, COMPLEX *, COMPLEX *, COMPLEX *); 
% main() 
% { 
% int n, k; 
% COMPLEX x[NUMPOINTS+1]; /* array of complex DFT data */ 
% COMPLEX A[NUMPOINTS]; /* array of complex A coefficients */ 
% COMPLEX B[NUMPOINTS]; /* array of complex B coefficients */ 
% COMPLEX IA[NUMPOINTS]; /* array of complex A* coefficients */ 
% COMPLEX IB[NUMPOINTS]; /* array of complex B* coefficients */ 
% COMPLEX G[2*NUMPOINTS]; /* array of complex DFT result */ 
% for(k=0; k<NUMPOINTS; k++) 
for k=0:(NUMPOINTS-1) 
    % { 
    % A[k].imag = (short)(16383.0*(-cos(2*pi/(double)(2*NUMPOINTS)*(double)k))); 
    % A[k].real = (short)(16383.0*(1.0 - sin(2*pi/(double)(2*NUMPOINTS)*(double)k))); 
    % B[k].imag = (short)(16383.0*(cos(2*pi/(double)(2*NUMPOINTS)*(double)k))); 
    % B[k].real = (short)(16383.0*(1.0 + sin(2*pi/(double)(2*NUMPOINTS)*(double)k))); 
    % IA[k].imag = -A[k].imag; 
    % IA[k].real = A[k].real; 
    % IB[k].imag = -B[k].imag; 
    % IB[k].real = B[k].real; 
    % } 
    A(k+1, 2) = floor(16383.0*(-cos(2*pi/(2*NUMPOINTS)*k))); 
    A(k+1, 1) = floor(16383.0*(1.0 - sin(2*pi/(2*NUMPOINTS)*k))); 
    B(k+1, 2) = floor(16383.0*(cos(2*pi/(2*NUMPOINTS)*k))); 
    B(k+1, 1) = floor(16383.0*(1.0 + sin(2*pi/(2*NUMPOINTS)*k))); 
    IA(k+1, 2) = -A(k+1, 2); 
    IA(k+1, 1) = A(k+1, 1); 
    IB(k+1, 2) = -B(k+1, 2); 
    IB(k+1, 1) = B(k+1, 1); 
end 
% /* Forward DFT */ 
% /* From the 2N point real sequence, g(n), for the N-point complex sequence, x(n) */ 
% for (n=0; n<NUMPOINTS; n++) 
% { 
for n=0:(NUMPOINTS-1) 
    % x[n].imag = g[2*n + 1]; /* x2(n) = g(2n + 1) */ 
    % x[n].real = g[2*n]; /* x1(n) = g(2n) */ 
    % } 
    x(n+1,2)=g(2*n + 1+1); 
    x(n+1,1)=g(2*n +1); 
end 
% /* Compute the DFT of x(n) to get X(k) -> X(k) = DFT{x(n)} */ 
% dft(NUMPOINTS, x); 
% void dft(int N, COMPLEX *X) 
% { 
% int n, k; 
% double arg; 
% int Xr[1024]; 
% int Xi[1024]; 
% short Wr, Wi; 
% for(k=0; k<N; k++) 
% { 
N=NUMPOINTS; 
for k=0:(N-1) 
    % Xr[k] = 0; 
    % Xi[k] = 0; 
    Xr(k+1)=0; 
    Xi(k+1)=0; 
    % for(n=0; n<N; n++) 
    % { 
    for n=0:(N-1) 
     % arg =(2*PI*k*n)/N; 
     % Wr = (short)((double)32767.0 * cos(arg)); 
     % Wi = (short)((double)32767.0 * sin(arg)); 
     % Xr[k] = Xr[k] + X[n].real * Wr + X[n].imag * Wi; 
     % Xi[k] = Xi[k] + X[n].imag * Wr – X[n].real * Wi; 
     arg = (2*pi*k*n)/N; 
     Wr = floor(32767*cos(arg)); 
     Wi = floor(32767*sin(arg)); 
     Xr(k+1) = Xr(k+1)+x(n+1,1)*Wr+x(n+1,2)*Wi; 
     Xi(k+1) = Xr(k+1)+x(n+1,2)*Wr-x(n+1,1)*Wi; 
     % } 
     % } 
    end 
end 
% for (k=0;k<N;k++) 
% { 
for k=0:(N-1) 
    % X[k].real = (short)(Xr[k]>>15); 
    % X[k].imag = (short)(Xi[k]>>15); 
    x(k+1,1)=floor(Xr(k+1)/pow2(15)); 
    x(k+1,2)=floor(Xi(k+1)/pow2(15)); 
    % } 
    % } 
end 
% /* Because of the periodicity property of the DFT, we know that X(N+k)=X(k). */ 
% x[NUMPOINTS].real = x[0].real; 
% x[NUMPOINTS].imag = x[0].imag; 
x(NUMPOINTS+1,1)=x(1,1); 
x(NUMPOINTS+1,2)=x(1,2); 
% /* The split function performs the additional computations required to get 
% G(k) from X(k). */ 
% split(NUMPOINTS, x, A, B, G); 
% void split(int N, COMPLEX *X, COMPLEX *A, COMPLEX *B, COMPLEX *G) 
% { 
% int k; 
% int Tr, Ti; 
% for (k=0; k<N; k++) 
% { 
for k=0:(NUMPOINTS-1) 
    % Tr = (int)X[k].real * (int)A[k].real – (int)X[k].imag * (int)A[k].imag + 
    % (int)X[N–k].real * (int)B[k].real + (int)X[N–k].imag * (int)B[k].imag; 
    Tr = x(k+1,1)*A(k+1,1)-x(k+1,2)*A(k+1,2)+x(NUMPOINTS-k+1,1)*B(k+1,1)+x(NUMPOINTS-k+1,2)*B(k+1,2); 
    % G[k].real = (short)(Tr>>15); 
    G(k+1,1)=floor(Tr/pow2(15)); 
    % Ti = (int)X[k].imag * (int)A[k].real + (int)X[k].real * (int)A[k].imag + 
    % (int)X[N–k].real * (int)B[k].imag – (int)X[N–k].imag * (int)B[k].real; 
    Ti = x(k+1,2)*A(k+1,1)+x(k+1,1)*A(k+1,2)+x(NUMPOINTS-k+1,1)*B(k+1,2)-x(NUMPOINTS-k+1,2)*B(k+1,1); 
    % G[k].imag = (short)(Ti>>15); 
    G(k+1,2)=floor(Ti/pow2(15)); 
    % } 
end 
% } 
% /* Use complex conjugate symmetry properties to get the rest of G(k) */ 
% G[NUMPOINTS].real = x[0].real - x[0].imag; 
% G[NUMPOINTS].imag = 0; 
% for (k=1; k<NUMPOINTS; k++) 
% { 
% G[2*NUMPOINTS-k].real = G[k].real; 
% G[2*NUMPOINTS-k].imag = -G[k].imag; 
% } 
G(NUMPOINTS+1,1) = x(1,1) - x(1,2); 
G(NUMPOINTS+1,2) = 0; 
for k=1:(NUMPOINTS-1) 
    G(2*NUMPOINTS-k+1,1) = G(k+1,1); 
    G(2*NUMPOINTS-k+1,2) = -G(k+1,2); 
end 
for k=1:(NUMDATA) 
    H(k)=sqrt(G(k,1)*G(k,1)+G(k,2)*G(k,2)); 
end 
end 

Еще одно:

function [ fr, fi ] = fix_fft(fr, fi) 
%UNTITLED Summary of this function goes here 
% Detailed explanation goes here 
N_WAVE = 1024; % full length of Sinewave[] 
LOG2_N_WAVE = 10;  % log2(N_WAVE) 
m = nextpow2(length(fr)); 

% void fix_fft(short fr[], short fi[], short m) 
% { 
% long int mr = 0, nn, i, j, l, k, istep, n, shift; 
mr=0; 
% short qr, qi, tr, ti, wr, wi; 
% 
% n = 1 << m; 
n = pow2(m); 
% nn = n - 1; 
nn = n-1; 
% 
% /* max FFT size = N_WAVE */ 
% //if (n > N_WAVE) return -1; 
% 
% /* decimation in time - re-order data */ 
% for (m=1; m<=nn; ++m) 
for m=1:nn 
    % { 
    %  l = n; 
    l=n; 

    %  do 
    %  { 
    %   l >>= 1; 
    %  } while (mr+l > nn); 
    not_done = true; 
    while(mr+l>nn || not_done) 
     l=floor(l/2); 
     not_done=false; 
    end 

    % 
    %  mr = (mr & (l-1)) + l; 
    mr = (mr & (l-1)) + l; 
    %  if (mr <= m) continue; 
    if (mr <= m) 
     continue 
    end 
    % 
    %  tr = fr[m]; 
    %  fr[m] = fr[mr]; 
    %  fr[mr] = tr; 
    %  ti = fi[m]; 
    %  fi[m] = fi[mr]; 
    %  fi[mr] = ti; 

    tr = fr(m+1); 
    fr(m+1) = fr(mr+1); 
    fr(mr+1) = tr; 
    ti = fi(m+1); 
    fi(m+1) = fi(mr+1); 
    fi(mr+1) = ti; 
    % } 
end 
% 
% l = 1; 
% k = LOG2_N_WAVE-1; 
l=1; 
k = LOG2_N_WAVE-1; 
% 
% while (l < n) 
% { 
while (l < n) 
    %  /* 
    %   fixed scaling, for proper normalization -- 
    %   there will be log2(n) passes, so this results 
    %   in an overall factor of 1/n, distributed to 
    %   maximize arithmetic accuracy. 
    % 
    %   It may not be obvious, but the shift will be 
    %   performed on each data point exactly once, 
    %   during this pass. 
    %  */ 
    % 
    %  // Variables for multiplication code 
    %  long int c; 
    %  short b; 
    % 
    %  istep = l << 1; 
    istep = l*2; 
    %  for (m=0; m<l; ++m) 
    %  { 
    for m=0:(l-1) 
     %   j = m << k; 
     %   /* 0 <= j < N_WAVE/2 */ 
     %   wr = Sinewave[j+N_WAVE/4]; 
     %   wi = -Sinewave[j]; 
     j = m*(pow2(k)); 
     wr = sin((j+N_WAVE/4)*2*pi/N_WAVE)*32768; 
     wi = -sin(j*2*pi/1024)*32768; 
     % 
     %   wr >>= 1; 
     %   wi >>= 1; 
     wr = floor(wr/2); 
     wi = floor(wi/2); 
     % 
     %   for (i=m; i<n; i+=istep) 
     %   { 
     i=m; 
     while(i<n) 
      %    j = i + l; 
      j = i+l; 
      % 
      %    // Here I unrolled the multiplications to prevent overhead 
      %    // for procedural calls (we don't need to be clever about 
      %    // the actual multiplications since the pic has an onboard 
      %    // 8x8 multiplier in the ALU): 
      % 
      %    // tr = FIX_MPY(wr,fr[j]) - FIX_MPY(wi,fi[j]); 
      %    c = ((long int)wr * (long int)fr[j]); 
      %    c = c >> 14; 
      %    b = c & 0x01; 
      %    tr = (c >> 1) + b; 
      c = wr * fr(j+1); 
      c = floor(c/pow2(14)); 
      b = c & 1; 
      tr = floor(c /2) + b; 
      % 
      %    c = ((long int)wi * (long int)fi[j]); 
      %    c = c >> 14; 
      %    b = c & 0x01; 
      %    tr = tr - ((c >> 1) + b); 
      c = wi * fi(j+1); 
      c = floor(c/pow2(14)); 
      b = c & 1; 
      tr = tr - (floor((c/2)) + b); 
      % 
      %    // ti = FIX_MPY(wr,fi[j]) + FIX_MPY(wi,fr[j]); 
      %    c = ((long int)wr * (long int)fi[j]); 
      %    c = c >> 14; 
      %    b = c & 0x01; 
      %    ti = (c >> 1) + b; 
      c = wr*fi(j+1); 
      c = floor(c/pow2(14)); 
      b = c & 1; 
      ti = floor((c /2)) + b; 
      % 
      %    c = ((long int)wi * (long int)fr[j]); 
      %    c = c >> 14; 
      %    b = c & 0x01; 
      %    ti = ti + ((c >> 1) + b); 
      c = wi * fr(j+1); 
      c = floor(c/pow2(14)); 
      b = c & 1; 
      ti = ti + (floor((c/2)) + b); 
      % 
      %    qr = fr[i]; 
      %    qi = fi[i]; 
      %    qr >>= 1; 
      %    qi >>= 1; 
      % 
      %    fr[j] = qr - tr; 
      %    fi[j] = qi - ti; 
      %    fr[i] = qr + tr; 
      %    fi[i] = qi + ti; 
      qr = fr(i+1); 
      qi = fi(i+1); 
      qr = floor(qr/2); 
      qi = floor(qi/2); 

      fr(j+1) = qr - tr; 
      fi(j+1) = qi - ti; 
      fr(i+1) = qr + tr; 
      fi(i+1) = qi + ti; 
      %   } 
      i = i+istep; 
     end 
     %  } 
    end 
    % 
    %  --k; 
    %  l = istep; 
    % } 
    k=k-1; 
    l=istep; 
end 
% } 
end 

Те в комментариях являются оригинальными, те не переведенный код

Тогда я моделируется с этой

function [ r ] = mfft(f) 
%MFFT Summary of this function goes here 
% Detailed explanation goes here 
Fs = 2048; 
T = 1/Fs; 
L = 2048; 
t = (0:L-1)*T; 
NFFT = 2^nextpow2(L); 
l = length(f); 
y = 0; 
for k=1:l 
    y = y + sin(2*pi*f(k)*t); 
end 
%sound(y, Fs); 
Y = fft(y,NFFT)/L; 
YY = fft_2(y)/L; 
[Y1 Y2] = fix_fft(y, zeros(1, L)); 
YYY = Y1+j()*Y2; 
f = Fs/2*linspace(0,1,NFFT/2+1); 
plot(f, 2*abs(Y(1:NFFT/2+1)), ':+b'); 
hold on 
plot(f, 2*abs(YY(1:NFFT/2+1)), '-or'); 
plot(f, abs(YYY(1:NFFT/2+1)), '--*g'); 
hold off 
r=0; 

end 

В принципе создать простую гармоническую волну с определенной частотой (скажем, 400Hz)

Ожидаемый выход БПФ должен быть шип на только 400, что соответствует встроенной функции БПФ, но другие коды не

Вот вывод графа

Синяя линия от встроено функции и что ожидается

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

зеленая линия является абсолютной беспорядок

Я попытался проверка программы несколько раз, но безрезультатно

я портирован их неправильно или каким-то образом я могу е Их?

Answer 1

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

В настоящее время FT не так медленно, как вы думаете. Если проект не должен трансформировать нечто вроде 10 000 точек данных менее чем за несколько миллисекунд. Кроме того, FT, по сравнению с FFT, прост и легко отлаживается.

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

Далее, используйте предварительно написанный пакет для БПФ. Сканирование в Интернете, я знаю, что есть пакеты, написанные на C, которые делают FFT.

В-третьих, если вам просто нужно написать свой собственный БПФ, тогда сделайте это. Но только если задачи (1) или (2) не соответствуют вашим требованиям. Будет сложно выполнить любые предварительно написанные пакеты FFT.

Вы узнаете многое на этом пути.