傅里叶matlab实现

参考：
https://blog.csdn.net/weixin_49082066/article/details/120926095
https://blog.csdn.net/weixin_42533243/article/details/117180030

1、文件

comp_minus.m

%u,v 复数
%u.a 实数， u.b虚数
function res = comp_minus(u, v)
e.a = u.a - v.a;
e.b = u.b - v.b;

res = e;

comp_plus.m

%u,v 复数
%u.a 实数， u.b虚数
function res = comp_plus(u, v)
e.a = u.a + v.a;
e.b = u.b + v.b;

res = e;

comp_times.m

%u,v 复数
%u.a 实数， u.b虚数
function res = comp_times(u, v)
e.a = u.av.a - u.bv.b;
e.b = u.av.b + u.bv.a;

res = e;

comp_impl.m

global N;
global U;

n = len/2;
if (len<2)
return; %Level 1
end

fftimpl(l,l+n-1, n);%even seg process
fftimpl(l+n, r, n);%odd seg process

for i = l : r
if i < (l + n)
wi = (i-l)(N/len);
in = i + n;
print_str = ['i= ',num2str(i), ’ wi = ’ ,num2str(wi), ’ in = ', num2str(in)];
disp(print_str);
cp = comp_plus(x(i + 1) ,comp_times(W(wi + 1) ,x(in + 1)));
U.a(i + 1) = cp.a;
U.b(i + 1) = cp.b;
else
in = i - n;
wi = (i-n-l)(N/len);
print_str = ['i= ',num2str(i), ’ wi = ',num2str(wi), ’ in = ', num2str(in)];
disp(print_str);
cp = comp_minus(x(in + 1),comp_times(x(i + 1),W(wi + 1)));
x(i + 1) = cp;
x(in + 1).a = U.a(in + 1);
x(in + 1).b = U.b(in + 1);
end
end

fftinit.m

global N;
N = 1024;

global fs;
fs = 2048;

global x;
%x(1).a = 0;
%x(1).b = 0;

global W;
%W = zeros(1 ,N);

global resolution;
resolution = fs/N;

rev.m

function res = rev(X,N)

len = log2(N);
u.a = zeros(1 ,N);
u.b = zeros(1 ,N);

%按位顺序反转序列
for i = 0 : N -1
idd = i;
ret = 0;
for j = 0 : len -1
%左移一位ret = ret << 1;
ret = bitshift(ret ,1);
% bit = idd & 1;
bit = bitand(idd ,1);
%idd = idd >> 1;
idd = bitshift(idd ,-1);
% ret = bit | ret;
ret = bitor(bit ,ret);
end
XX = ['i = ’ ,num2str(i) , ’ ret = ’ ,num2str(ret)];
%disp(XX);
u.a(i + 1) = X(ret + 1).a;
u.b(i + 1) = X(ret + 1).b;
end

for i = 1 : N
X(i).a = u.a(i);
X(i).b = u.b(i);
end

res = X;

Wn.m

%float a , float b
%pi 3.1415926
function ret = Wn(A ,B ,pi)
u.a = cos((2*pi/A)B);
u.b = -sin((2pi/A)*B);

ret = u;

测试验证testfftimpl.m

global W;
global U;
U.a = zeros(1 ,N);
U.b = zeros(1 ,N);

pi = 3.1415926;

for i = 1 : N
t=i/fs; %time resolution
%Generate the original signal sequence and tranfer it to Complex number field
x(i).a = 114sin(2pi50t) + 514sin(2pi152t) + 1919sin(2pi536t) + 810sin(2pi996t);
x(i).b=0;
%u.a = cos((2pi/N)i);
%u.b = -sin((2pi/N)i);
cp = Wn(N,i ,pi);
W(i).a = cp.a;%e^-j((2pii)/N)
W(i).b = cp.b;
end
x_list = cat(1 ,x.a);
figure(‘name’, ‘原始信号’); grid on; hold on;
plot(x_list);

x = rev(x,N);
fftimpl(0 ,N - 1 ,N);

P2 = zeros(1,N);
for i = 1 : N
P2(i) = sqrt(x(i).ax(i).a + x(i).bx(i).b)/N;
end
P1 = P2(1:N/2 + 1);
P1(2:end-1) = 2*P1(2:end-1);

Fs = 2048;
figure(‘name’, ‘傅里叶变换-代码实现’); grid on; hold on;
f = Fs*(0:(N/2))/N; %1751数组1000(0:750)/1500 [0 0.66 1.33 …
plot(f,P1)

Y = fft(x_list);
oP2 = abs(Y/N);
oP1 = oP2(1:N/2 + 1);
oP1(2:end-1) = 2*oP1(2:end-1);

figure(‘name’, ‘傅里叶变换-fft’); grid on; hold on;
f = Fs*(0:(N/2))/N; %1751数组1000(0:750)/1500 [0 0.66 1.33 …
plot(f,oP1)

原理：

所以：


根据这几个特性，就可以将一个长的DFT运算分解为若干短序列的DFT运算的组合，从而减少运算量




那么，通过这样的分解，每一个N/2点DFT只需(N^2)/4次复数相乘计算，明显的节省了运算量。

那么以此类推，继续将已经得出的X1(k)和X2(k)按奇偶继续分解，还是上边一样的方法。那么就得出了百度上都可以找到的一大堆的这个图片了。