基于MATLAB去理解掌握傅里叶级数和傅里叶变换

1.周期信号的傅里叶级数

f(t)=f(t+T) F0=1/T为基波频率
满足狄利赫里条件则周期信号可以展开为三角函数的线性组合
（1） 在一个周期内，函数f（t）为连续或只含有有限个第一类间断点；
（2） 在一个周期内，函数f（t）的极值点为有限个；
（3） 在一个周期内，函数f（t）是绝对可积的；

利用欧拉公式将其转换为复指数形式的傅里叶级数

三角函数标准形式中cn是第n次谐波分量的振幅，但在指数形式中，Fn要与相对应的F–n合并，构成第n次谐波分量的振幅和相位。是cn幅度一半

2.周期信号的频谱

一个周期信号与另一个周期信号的区别，在时域中表现为波形不同；而在频域中则表现为Fn不同，即振幅和相位的不同。因而振幅和相位是在频域中研究信号f（t）的关键。

前面已述，周期信号的复振幅Fn一般为nω0的复函数，因而描述其特点的频谱图一般有两个：一个称为振幅频谱，简称幅度谱，它是以ω为横坐标、振幅为纵坐标所画的谱线图；另一个称为相位频谱，简称相位谱，它是以ω为横坐标、相位为纵坐标所画的谱线图。

在信号的复振幅Fn为ω的实函数的特殊情况下，其复振幅与变量（ω）的关系也可以用一个图绘出。

信号的时域波形与频谱都是实际存在的，例如我们可以通过示波器来观察信号的时域波形，通过频谱分析仪观察信号的频谱。声波有频谱，图像也有频谱，频谱与时域波形一样具有实际意义。

试画出周期信号f（t）=–1+2sin（0.2πt）–3cosπt的幅度频谱。
解：MATLAB程序如下：

% %周期信号的幅度频谱 
t=0:0.1:50;
f=-1+2*sin(0.2*pi*t)-3*cos(pi*t);
F=fftshift(fft(f))/501;
subplot(211)
plot(t,f);
axis([-1 50 -6.1 4.1]);
subplot(212) 
w=-5*pi:0.02*pi:5*pi; 
plot(w,abs(F))%注意w和F的点数要一样 直接看点点数求出来，以下面的一种方法为标准方法
axis([-5 5 -0.1 3]);
 
figure
N=5000;T=0.1;n=1:N;  %5000个数据点 时间间隔 0.1s
D=2*pi/(N*T);  %将时间域转换成频域  频域间隔0.004pi 因为fft横轴代表的是点数，所以要乘以D来转换成角频率
f=-1+2*sin(0.2*pi*n*T)-3*cos(pi*n*T);
F=fftshift(fft(f))/N;     %为什么要除以N，因为DFT的频谱除以信号长度才是真实的周期信号傅里叶级数值
%又因为是双边谱，所以0.2pi是1 pi是1.5
k=floor(-(N-1)/2:N/2);    %将正域的5000个点转移至正负域，k都是整数
subplot(311)
plot(n*T,f);
axis([-1 50 -6.1 4.1]);
line([-1,50],[0,0]);
line([0,0],[-6.1,4.1]);
subplot(312) 
plot(k*D,abs(F))%注意w和F的点数要一样
axis([-6 6 -0.1 3]);
subplot(313) 
plot(k*D,angle(F))%注意w和F的点数要一样
axis([-1 1 -4 4]);
 

构造基频为1Hz的方波信号分析其傅里叶级数
注意这个方波的形式，为+1 -1 τ=0.5 T=1； 所以不含直流分量

N=5000;Ts=0.001;n=1:N;   %时间间隔 0.1s 选好时间向量
f=square(2*pi*n*Ts);  %产生方波
F=fftshift(fft(f))/N;
k=floor(-(N-1)/2:N/2);
subplot(211)
plot(n*Ts,f)
axis([0,2,-2,2]);
subplot(212)
f0=1/(N*Ts) %单位Hz
plot(k*f0,abs(F))
axis([-20,20,0,1]);

可以看出无直流分量，在基频出幅度为4*1/π（单边普）2/π（双边普）

区别fourier变换和fft

clear
N=5000;Ts=0.1;n=1:N;
f0=1/(N*Ts);%频率分辨率 单位hz
fs=1/Ts;
k=floor(-(N-1)/2:N/2);
w=k*f0;
t=n*Ts;
f=cos(pi*t);
F=T*f*exp(-j*t'*w);  %fourier变换，其横轴就是角频率
subplot(211)
plot(w,abs(F));   %pi
subplot(2,1,2);
F1=fftshift(fft(f))/N;
plot(w,abs(F1));%fft横轴是点数k，然后转换成k*f0,所以是在0.5hz处，若使f0*2pi就转换成加频率了

这时的fourier的横轴是1/NTs=1/500 为频率间隔 在pi处有频谱泄漏，因为Ts=0.1 不是整周期取，将频率间隔变成2pi/50，在pi处不泄露到其它地方
而fft满足抽样定理，在那个点上是不会泄露的就是真实的幅度

与点数有关 n所表示的频率为：Fn=(n-1)*Fs/N。由上面的公式可以看出，Fn所能分辨到频率为为Fs/N。如果采样频率Fs为1024Hz，采样点数为1024点，则可以分辨到1Hz。

clear
N=5000;T=0.1;n=1:N;
f0=2*pi/(N*T);%频率分辨率 
fs=1/T;
k=floor(-(N-1)/2:N/2);
w=k*f0;
t=n*T;
f=cos(pi*t);
F=f*exp(-j*t'*w)/N;  %fourier变换，其横轴就是角频率
subplot(211)
plot(w,abs(F));   %pi
subplot(2,1,2);
F1=fftshift(fft(f))/N;
plot(w,abs(F1));%fft横轴是点数k，然后转换成k*f0,所以是在0.5hz处，若使f0*2pi就转换成加频率了

非周期信号的傅里叶变换

由于基波频率趋于无穷小量，因此非周期信号f（t）包含了所有的频率分量，即非周期信号的频谱为连续谱。由于各频率成分的振幅趋于无穷小，因此非周期信号f（t）的频谱只能用密度函数|F（jω）|来表述各分量的相对大小。
注意要是fourier在频域上出现dirac函数，无法在图上显示 无穷大

双边指数函数f（t）=e–a|t|，其中–∞＜t＜∞，a＞0
或f（t）=eatε（–t）+e–atε（t）
利用以上单边指数函数傅里叶变换的结果，有

clear
syms t v;
F=fourier(exp(-2*abs(t))) %傅里叶变换就是对非周期信号求得其F(w)
subplot(2,1,1);
ezplot(exp(-2*abs(t)));
subplot(212);
ezplot(F);

画出尺度变换后的门函数的频谱图，f1（t）=ε（t+1/2）–ε（t–1/2）与f2（t）=ε（t+2）–ε（t–2），写出对应的傅里叶变换的MATLAB程序及程序运行结果。
常用此数值法求得非周期信号的傅里叶变换的频谱

clear;
 T=0.02;t=-10:T:10;N=200;
 W=4*pi;k=-N:N;w=k*W/N; 
 f1=(stepfun(t,-0.5)-stepfun(t,0.5));  %f1(t)??  
 f2=(stepfun(t,-2)-stepfun(t,2));  %f1(t)?
 F=T*f1*exp(-j*t'*w);            %f(t)的傅里叶变换?? 数组积分求面积傅里叶变换
 F1=abs(F);P1=angle(F); 
 FF=T*f2*exp(-j*t'*w);  
 F2=abs(FF);P2=angle(FF); 
 
 subplot(3,2,1);plot(t,f1);  
 axis([-3,3,-1.2,1.2]);ylabel('f1(t)');
 xlabel('t');title('f1(t)');grid;
 
 subplot(3,2,2);plot(t,f2);  
 axis([-3,3,-1.2,1.2]);ylabel('f2(t)');
 xlabel('t');title('f2(t)');grid;
 
 subplot(3,2,3);plot(w,F1);
 axis([-3*pi,3*pi,-0.01,2.1]);
 grid;ylabel('振幅');
 
 subplot(3,2,4);plot(w,F2);
 axis([-3*pi,3*pi,-0.01,4.1]);
 grid;ylabel('振幅');
 
 subplot(3,2,5);plot(w,P1*180/pi);
 grid;axis([-3*pi,3*pi,-180,180]);
 xlabel('w');ylabel('相位(度)')
 
 subplot(3,2,6);plot(w,P2*180/pi);
 grid;axis([-3*pi,3*pi,-180,180]);
 xlabel('w');ylabel('相位(度)')

写在后面

这些基本上包括了傅里叶级数和傅里叶变换在matlab中的体现，更详细更多的例子详见word文档，后续会更新dsp的内容噢！