我希望在微控制器上实现 FFT 算法，所以我想在实际使用之前模拟代码
我有 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）的简单正弦波
FFT 的预期输出应该只是 400 处的尖峰，内置 FFT 函数同意，但其他代码不同意
这是输出图
蓝线来自内置函数，是预期的
红线是上面的代码，看起来不错，除了在其他地方有更高幅度的尖峰
绿线绝对是一团糟
我尝试检查程序几次但无济于事
我是不是把它们移植错了，还是我可以修复它们？