实验三连续周期性时间信号的傅里叶级数
一、实验目的：
1. 进一步掌握MATLAB子函数的表示方法
2. 深刻理解傅里叶级数的信号分解理论及收敛性问题
3. 理解周期性信号的频谱特点。
二、实验原理
傅里叶级数
设有连续时间周期信号，它的周期为T，角频率，且满足狄里赫利条件，则该周期信号可以展开成傅里叶级数，即可表示为一系列不同频率的正弦或复指数信号之和。
傅里叶级数有三角形式和指数形式两种。
1. 三角形式的傅里叶级数：
式中系数，称为傅里叶系数，可由下式求得：
[
2. 指数形式的傅里叶级数：
式中系数称为傅里叶复系数，可由下式求得：
周期信号频谱具有三个特点：
（1）离散性，即谱线是离散的；
（2）谐波性，即谱线只出现在基波频率的整数倍上；
（3）收敛性，即谐波的幅度随谐波次数的增高而减小。
周期信号的MATLAB表示
周期信号的傅里叶分解用Matlab进行计算时，本质上是对信号进行数值积分运算。
在Matlab中有多种进行数值积分运算的方法，我们采用quadl函数，它有两种其调用形式。
(1) y＝quadl(‘func’, a, b)。
其中func是一个字符串，表示被积函数的.m文件名（函数名）；a、b分别表示定积分的下限和上限。
(2) y＝quadl(@myfun, a, b)。
其中“@”符号表示取函数的句柄，myfun表示所定义函数的文件名。
例：
用MATLAB计算脉冲宽度T1 = 2；周期T = 4的周期性脉冲信号的复傅里叶级数，分别画出N = -2:2， -10:10， -50:50， -200:200的傅里叶级数展开及合成，观察吉普斯效应。
画出T = 4， T =8下的双边谱
A.首先创建一个子函数singRect（t, T1），表示单个脉冲信号，时间为t，宽度为T1。
function y = singRect(t, T1)
y = (abs(t) <= T1);
end
B.创建傅里叶积分的被积子函数
function y = rectExp(t, k, w)
y = (abs(t) <= 1) .* exp(-1j*k*w*t);
end
C.创建子函数用于傅里叶级数计算及合成
function [x, ak] = fourierSeries(N, t)
T1 = 1;
T = 4; w = 2 * pi/T;
ak = zeros(1, 2 * N + 1);
for i = 1:2*N+1 %傅里叶分解，计算傅里叶系数ak
ak(i) = quadl(@(t)fsInt(t, i - N - 1, w, T1), -2, 2)/T;
end;
x = 0;
for i = 1:2*N + 1 %傅里叶级数合成
x = x + ak(i) * exp(1j*(i - N - 1)*w*t);
end
end
D.创建main函数，计算不同N下的傅里叶级数及合成。
T1 = 1;
T = 4;
t = -T/2:0.001:T/2;
figure,
subplot 221,
N = 2;
[x, ak] = fourierSeries(N, t);
plot(t, singRect(t, T1), 'k');
hold on ;
plot(t, real(x), 'r'); title('N = 2'); xlim([-2 2]);
subplot 222, N = 10;
[x, ak] = fourierSeries(N, t); plot(t, singRect(t, T1), 'k'); hold on ;
plot(t, real(x), 'r'); title('N = 10'); xlim([-2 2]);
subplot 223, N = 50;
[x, ak] = fourierSeries(N, t); plot(t, singRect(t, T1), 'k'); hold on ;
plot(t, real(x), 'r'); title('N = 50'); xlim([-2 2]);
subplot 224, N = 200;
[x, ak] = fourierSeries(N, t); plot(t, singRect(t, T1), 'k'); hold on ;
plot(t, real(x), 'r'); title('N = 300'); xlim([-2 2]);
运行结果如下图所示，可以看到明显的Gibbs 效应，随着N 的增加，合成越接近于原信号，但Gibbs 效应仍然存在。
N = 2
N = 10
N = 50
N = 300
T1 = 1, T = 8
T1 = 1, T = 4
上图是不同T下的频谱，可以看到T越大，谱线越密。
四、作业
计算连续周期性三角波的傅里叶级数，设周期为 T = 5s，画出N = 5；N = 10；N = 100；N = 300时的傅里叶级数合成图。
并画出T = 5,10,20的傅里叶系数a k。
function y=FsIntegralFunc(t,k,T,width)
w=2*pi/T
y=tripuls(t,width,0).*exp(-1j*k*w*t);
end
function ak=fourierSeries(t,k,T,width)
ak=quadl(@(t)FsIntegralFunc(t,k,T,width),-T/2,T/2)/T;
end
N=input('please Input N= ');
k=-N:N
width=4;
T=5;
t=-T/2:0.01:T/2;
w=2*pi/T;
ak=zeros(1,length(k));
x=0;
N=10
figure;
for m=-N:N;
ak(m+N+1)=fourierSeries(t,m,T,width);
x= x + ak(m+N+1) * exp(1j*m*w*t);
end
stem(k,real(ak),'k');
title('傅里叶系数ak')
xlabel('T=5');
figure('name','信号的合成');
plot(t,tripuls(t,width,0));
hold on;
plot(t,x,'r');
title(‘N= ’)
T=5
傅里叶系数ak
T=10
傅里叶系数ak
T=20。