快速傅里叶变换
%快速傅里叶变换
function ret_val = myfft(Vector)
%因为输入的数据可能不是2的整数次幂,变换使得计算更加方便
[m,n]=size(Vector);%输入信号矩阵大小
num=ceil(log2(n));%向上取整
N=2^num;
vector=zeros(m,N);%申请足够大小矩阵
vector(:,1:n)=Vector(:😅;%将变换后的信号输入,不足补零
%变址运算
% temp=zeros(1,N);
% for p=1:m
% for q=0:N-1
% t=bin2dec(fliplr(dec2bin(q,num)));%变址运算
% temp(1,q+1)=vector(p,t+1);
% end
% vector(p,:)=temp(1,:);
% end
%变址运算
for line=1:m%循环行数
j1 = 0;
for i = 1 : N%循环个数
if i < j1 + 1
tmp = vector(line,j1 + 1);
vector(line,j1 + 1) = vector(line,i);
vector(line,i) =tmp;
end
k = N / 2;
while k <= j1
j1 = j1 - k;
k = k / 2;
end
j1 = j1 + k;
end
end
for lines=1:m%循环行数
for L=0:num-1%L表示运算等级或者层数
dis=2^L;%dis表示奇偶组之间的距离
for id=1:2^(num-L-1) %循环当前层数组数
%进行同址运算
for idx=1:dis%循环组内个数
x1=(id-1)2dis+idx;%求得奇数数组的索引值
x2=(id-1)2dis+dis+idx;%对应偶数数组的索引值
temp1=vector(lines,x1)+vector(lines,x2)W(L,(x1-1));%中间变量保存相应奇偶数组数据
temp2=vector(lines,x1)-vector(lines,x2)W(L,(x1-1));
vector(lines,x1)=temp1;%所谓存入之前地址
vector(lines,x2)=temp2;
end
end
end
end
ret_val =vector;
function val=W(L,x)%旋转因子当层数为L,索引值为x
val=exp(-1j2pi*x/2^(L+1));
end
end
例题
Fs = 100; % Sampling frequency
t = -0.5:1/Fs:0.5; % Time vector
L = length(t); % Signal length
X = 1/(4sqrt(2pi0.01))(exp(-t.^2/(2*0.01)));
plot(t,X)
title(‘Gaussian Pulse in Time Domain’)
xlabel(‘Time (t)’)
ylabel(‘X(t)’)
figure
n = 2^nextpow2(L);
Y = fft(X,n);
f = Fs*(0:(n/2))/n;
P = abs(Y/n).^2;
plot(f,P(1:n/2+1))
title(‘Gaussian Pulse in Frequency Domain’)
xlabel(‘Frequency (f)’)
ylabel(‘|P(f)|^2’)
傅里叶级数逼近
function [A0,A,B]=myfzz(func,T,n)
syms t;
func=subs(sym(func),sym(‘x’),sym(‘t’));
A0=int(sym(func),t,-T/2,T/2)/T;
for(k=1:n)
A(k)=int(funccos(2pikt/T),t,-T/2,T/2)*2/T;
A(k)=vpa(A(k),4);
B(k)= int(funcsin(2pikt/T),t,-T/2,T/2)*2/T;
B(k)= vpa(B(k),4);
end
例题