1、定义
T=1;% 周期0-1
N=20;% 最大谐波
k=-N:N;% -20:20——-20表示a_(20)e^(j*20*w0*t)
N1=length(k);%N1=41
% ceil(N1/2)=21%即21是N1的中位数
t=linspace(0,T,100);
Nt=length(t);
tt=linspace(-4*T,4*T,1024);
Ntt=length(tt);
2、原函数
w0=2*pi/T;
% 0-T
xt=(t>=T/4).* 1.0;
fig10=figure(10);
fig10.Name='原函数';
plot(t,xt);
axis([t(1) t(Nt) min(xt)-0.05 max(xt)+0.05 ]);
3、原函数的周期延拓
N2=length(tt);
xx=zeros(1,N2);%size(a)=1*41，用于记录傅里叶级数的系数
for n=1:N2
for numT=-4:3
if (T/4+numT*T)<tt(n) && tt(n)<(T+numT*T)
xx(n)=1;
break
end
end
end
fig20=figure(20);
fig20.Name='原函数的周期延拓';
plot(tt,xx);
axis([tt(1) tt(Ntt) min(xx)-0.05 max(xx)+0.05 ]);
grid on 
4、使用integral进行傅里叶变换
a=zeros(1,N1);%size(a)=1*41，用于记录傅里叶级数的系数
for n=1:N1
% a(n)=integral(@(x) signal0(x,k(n),T),0,T); %1
a(n)=integral(@(x) signal0(x,k(n),T),0,T,'RelTol',0,'AbsTol',1e-12);
% %解决《警告: 最小步长已快达到 x = 0.25。可能具有奇异性，或者容差可能对于此问题太小。 》
% RelTol - 相对误差容限
% AbsTol - 绝对误差容限
end
A=abs(a); %计算傅里叶变换的幅值
fig30=figure(30);
fig30.Name='傅里叶分解图';
subplot(3,1,1);stem(k,A);ylabel("a(k)幅值");
axis([k(1) k(N1) min(A)-0.05 max(A)+0.05]);% axis([ x起点 x终点 y起点 y终点]);
subplot(3,1,2); stem(k,real(a));ylabel('a(k)实部'); % real(F)
axis([k(1) k(N1) min(real(a))-0.05 max(real(a))+0.05 ]);
subplot(3,1,3);stem(k,imag(a));ylabel('a(k)虚部'); % imag(F)
xlabel('k');
axis([k(1) k(N1) min(imag(a))-0.05 max(imag(a))+0.05 ]);
5、谐波分量
%基波分量
fundamental_component=a(ceil(N1/2))*exp(1i*0*w0*tt);
%一次谐波分量
% first_harmonic_component=a(ceil(N1/2)-1)*exp(-1i*1*w0*tt)+a(ceil(N1/2)+1)*exp(1i*1*w0*tt);
har=1;
first_harmonic_component=a(ceil(N1/2)-har)*exp(-1i*har*w0*tt)+a(ceil(N1/2)+har)*exp(1i*har*w0*tt);
%二次谐波分量
% second_harmonic_component=a(ceil(N/2)-2)*exp(-1i*2*w0*tt)+a(ceil(N/2)+2)*exp(1i*2*w0*tt);
har=2;
second_harmonic_component=a(ceil(N1/2)-har)*exp(-1i*har*w0*tt)+a(ceil(N1/2)+har)*exp(1i*har*w0*tt);
%三次谐波分量
har=3;
third_harmonic_component=a(ceil(N1/2)-har)*exp(-1i*har*w0*tt)+a(ceil(N1/2)+har)*exp(1i*har*w0*tt);
fig40=figure(40);
fig40.Name='傅里叶分解图';
subplot(3,1,1);
plot(tt,fundamental_component);
ylabel('基波分量');
subplot(3,1,2);
plot(tt,fundamental_component+first_harmonic_component+second_harmonic_component);
ylabel(sprintf('%d%s',2,'次谐波分量合成图'));
subplot(3,1,3);
plot(tt,fundamental_component+first_harmonic_component+second_harmonic_component+third_harmonic_component)
ylabel(sprintf('%d%s',3,'次谐波分量合成图'));
6、任意次谐波的傅里叶级数合成
xiebo=8;% 计算前xiebo次谐波分量的和
% tt=linspace(-4*T,4*T,1024);
fx=a(ceil(N1/2))*exp(1i*0*w0*tt);
for n=1:xiebo
fx=fx+a(ceil(N1/2)-n)*exp(-1i*n*w0*tt)+a(ceil(N1/2)+n)*exp(1i*n*w0*tt);
% n=1——第20次谐波分量
% a(1)*exp(1i*(1-20-1)*w0*tt)=a(1)*exp(1i*(-20)*w0*tt)
% n=N1=41——第20次谐波分量
% a(41)*exp(1i*(41-20-1)*w0*tt)=a(41)*exp(1i*(20)*w0*tt)
end
fig41=figure(41);
fig41.Name=sprintf('%d%s',xiebo,'次谐波分量合成图');
% f3=figure('Name','傅里叶合成图','IntegerHandle','on','Number',100)
plot(tt,fx)
grid on
7、20次谐波的傅里叶级数合成
% tt=linspace(-4*T,4*T,1024);
fx=0;
for n=1:N1
fx=fx+a(n)*exp(1i*(n-N-1)*w0*tt);
% n=1——第20次谐波分量
% a(1)*exp(1i*(1-20-1)*w0*tt)=a(1)*exp(1i*(-20)*w0*tt)
% n=N1=41——第20次谐波分量
% a(41)*exp(1i*(41-20-1)*w0*tt)=a(41)*exp(1i*(20)*w0*tt)
end
fig42=figure(42);
fig42.Name='傅里叶合成图';
% f3=figure('Name','傅里叶合成图','IntegerHandle','on','Number',100)
plot(tt,fx)
grid on
8、局部函数——signal0
%程序signal0.m:定义求傅里叶级数的系数的被积函数: y=(1/T)*x(t).*exp(-j*(2*pi/T)*k*t)
function y=signal0(t,k,T) %T为信号周期，t为时间变量,k为k次谐波
x=(t>=T/4).* 1.0;
% 定义占空比为25%的非对称周期方波for one period=T
% 对称周期三角形可定义:x=(abs(t)<=T/4). * (1-abs(t))
% 占空比为50%的对称周期方波可定义:x=(abs(t)<=T/4). *1.0
y=(1/T).*x.*exp(-1i*k*(2*pi/T)*t);
end傅里叶分析和滤波- MATLAB & Simulink- MathWorks 中国
