标签：
1、掌握傅立叶级数（FS），学会分析连续时间周期信号的频谱分析及MATLAB实现；
2、掌握傅立叶变换（FT），了解傅立叶变换的性质以及MATLAB实现。
F1=int(f,v,a,b) — 对f表达式的v变量在（a,b)区间求定积分 F2=subs(s,OLD,NEW)－用新变量NEW代替S中的指定变量OLD。 F3=vpa(x,n) ： 显示可变精度计算；x为符号变量，n表示要精确计算的位数。
代码抄袭如下：
%ex_1 %求系数 clear all; syms t x n t0; T=10; % 信号周期 tao_2=0.5; %脉冲宽度 Nf=7; % 分解的最高级数 Nn=6; % 有效位数 x=heaviside(t+t0)-heaviside(t-t0); % 注意：Ver 2011b x=subs(x,t0,tao_2) A0=int(x,t,-tao_2,T-tao_2)/T % a0 As=int(x*2*cos(2*pi*n*t/T)/T,t,-tao_2,T-tao_2) % an Bs=int(x*2*sin(2*pi*n*t/T)/T,t,-tao_2,T-tao_2) % bn Fn=(As-j*Bs)/2 % fn A(1)=double(vpa(A0,Nn)); % 直流分量 for k=1:Nf A(k+1)=double(vpa(subs(As,n,k),Nn)); B(k+1)=double(vpa(subs(Bs,n,k),Nn)); end A B %求各次谐波：：开始采用数值处理方法 t1=-T/2:0.01:T/2; f0=A(1); %直流分量 f1=A(2).*cos(2*pi*1*t1/T)+B(2).*sin(2*pi*1*t1/T); ; % 基波 f2=A(3).*cos(2*pi*2*t1/T)+B(3).*sin(2*pi*2*t1/T); ; % 2次谐波 f3=A(4).*cos(2*pi*3*t1/T)+B(4).*sin(2*pi*3*t1/T); % 3次谐波 f4=A(5).*cos(2*pi*4*t1/T)+B(5).*sin(2*pi*4*t1/T); ; % 4次谐波 f5=A(6).*cos(2*pi*5*t1/T)+B(6).*sin(2*pi*5*t1/T); % 5次谐波 f6=A(7).*cos(2*pi*6*t1/T)+B(7).*sin(2*pi*6*t1/T); % 6次谐波 f7=f1+f2+f0; % 基波+2次谐波 f8=f7+f3+f0; % 基波+2次谐波+3次谐波 f9=f8+f4+f6+f0; % 基波+2次谐波+3次谐波+4次谐波+6次谐波 %画出谐波图形 y=subs(x,t,t1); %调用连续时间函数-周期矩形脉冲 subplot(2,2,1),plot(t1,f1),hold on;plot(t1,y,‘r:‘);title(‘周期矩形波的形成—基波‘),axis([-2.5,2.5,-0.5,1.1]) subplot(2,2,2),plot(t1,f7),hold on;plot(t1,y,‘r:‘);title(‘周期矩形波的形成—基波+2次谐波‘),axis([-2.5,2.5,-0.5,1.1]) subplot(2,2,3),plot(t1,f8),hold on;plot(t1,y,‘r:‘);title(‘基波+2次谐波+3次谐波‘),axis([-2.5,2.5,-0.5,1.1]) subplot(2,2,4),plot(t1,f9),hold on;plot(t1,y,‘r:‘);title(‘基波+2次谐波+3次谐波+4次谐波+6次谐波‘),axis([-2.5,2.5,-0.5,1.1]) %求频谱---------------这边给出双边谱的正轴部分，仅供参考 Fn=(As-j*Bs)/2 Nf=60; A(1)=double(vpa(A0,Nn)); for k=1:Nf A(k+1)=double(vpa(subs(As,n,k),Nn)); B(k+1)=double(vpa(subs(Bs,n,k),Nn)); end Fn1(1) = A(1); Fn1(2:Nf+1)=(A(2:Nf+1)-j*B(2:Nf+1))/2; absFn1 = abs(Fn1); angleFn1 = angle(Fn1); figure; subplot(211);stem([0:Nf],absFn1,‘r.‘);title(‘双边幅度谱的正半轴，注意与p200图比较幅度值‘) subplot(212);stem([0:Nf],angleFn1,‘r.‘);title(‘双边相位谱的正半轴‘)
注意：课件里的基波定义有误，不应该含有直流分量
% 例9_1 % 观察周期方波信号的分解与合成 % m：傅里叶级数展开的项数 clc;clear all;close all; display(‘Please input the value of m (傅里叶级数展开的项数)‘); % 在命令窗口显示提示信息 m = input(‘m = ‘); % 键盘输入傅里叶级数展开的项数 t = -2*pi:0.01:2*pi; % 时域波形的时间范围-2π~2π，采样间隔0.01 n = round(length(t)/4); % 根据周期方波信号的周期，计算1/2周期的数据点数 f = [ones(n,1);-1*ones(n,1);ones(n,1);-1*ones(n+1,1)]; %构造周期方波信号 y = zeros(m+1,max(size(t))); y(m+1,:) = f‘; figure(1); plot(t/pi,y(m+1,:)); %绘制方波信号 grid; %在图形中加入栅格 axis([-2 2 -1.5 1.5]); %指定图形显示的横坐标范围和纵坐标范围 title(‘周期方波‘); %给显示的图形加上标题 xlabel(‘单位pi‘,‘Fontsize‘, 8); %显示横坐标单位 x = zeros(size(t)); kk = ‘1‘; for k=1:2:2*m-1 %循环显示谐波叠加图形 pause; x = x+sin(k*t)/k; y((k+1)/2,:) = 4/pi*x; %计算各次谐波叠加和 plot(t/pi,y(m+1,:)); hold on; plot(t/pi,y((k+1)/2,:)); %绘制谐波叠加信号 hold off; grid; axis([-2 2 -1.5 1.5]); title(strcat(‘第‘,kk,‘次谐波叠加‘)); xlabel(‘单位pi‘,‘Fontsize‘, 8); kk = strcat(kk,‘、‘,num2str(k+2)); end pause; plot(t/pi,y(1:m+1,:)); grid; axis([-2 2 -1.5 1.5]); title(‘各次谐波叠加波形‘); xlabel(‘单位pi‘,‘Fontsize‘, 8); % End
%ex_3求频谱 clear all;syms t x n t0; T=5;tao=0.5;Nf=60;Nn=6; x=heaviside(t+t0)-heaviside(t-t0) x=subs(x,t0,tao) A0=int(x,t,-tao,T-tao)/T As=int(x*2*cos(2*pi*n*t/T)/T,t,-tao,T-tao) Bs=int(x*2*sin(2*pi*n*t/T)/T,t,-tao,T-tao) Fn=(As-j*Bs)/2 A(1)=double(vpa(A0,Nn)); for k=1:Nf A(k+1)=double(vpa(subs(As,n,k),Nn)); B(k+1)=double(vpa(subs(Bs,n,k),Nn)); end t1=-2.5:0.003:2.5; y=subs(x,t,t1); subplot(3,1,1),plot(t1,y),title(‘矩形脉冲‘),axis([-2.5,2.5,-0.1,1.1]) %单边谱 Fs(1)=A(1); Fs(2:Nf+1)=abs(A(2:Nf+1)-j.*B(2:Nf+1)); Ns=0:Nf; subplot(3,1,2),stem(Ns,Fs),title(‘连续时间函数的单边谱‘) %双边谱 N=Nf*2*pi/T; K=-N:2*pi/T:N; Fs(2:Nf+1)=abs(A(2:Nf+1)-j.*B(2:Nf+1))/2; Fd=[fliplr(Fs),Fs(2:Nf+1)]; subplot(3,1,3),stem(K,Fd),title(‘连续时间函数的双边谱‘)
% 周期方波信号频谱分析 function CTFS_SQ % 绘制并观察周期方波信号频谱 % Nf：傅里叶级数展开的项数 % an：各次谐波余弦项系数 display(‘Please input the value of Nf ‘); Nf = input(‘Nf = ‘); an = zeros(Nf+1,1); cn(1) = 0; for i = 1:Nf an(i+1) = 2/(i*pi)*sin(i*pi/2); %计算系数an cn(i+1) = abs(2/(i*pi)*sin(i*pi/2)); %计算幅度谱 end t = -5:0.01:5; x = square(pi*(t+0.5)); %用square函数求方波信号 subplot(211); plot(t,x); %绘制方波信号 axis([-5 5 -1.5 1.5]); title(‘周期方波信号‘,‘Fontsize‘,8); xlabel(‘t (单位:s)‘, ‘Fontsize‘,8); subplot(212); k = 0:Nf; stem(k,cn); hold on; plot(k,cn); title(‘幅度频谱‘,‘Fontsize‘,8); xlabel(‘谐波次数‘, ‘Fontsize‘,8); % End
% 周期锯齿脉冲信号频谱分析 function CTFS_JC % 绘制并观察周期锯齿脉冲信号频谱特性 % Nf：谐波的次数 % bn：第1,2,3,...次谐波正弦项系数 display(‘Please input the value of Nf ‘); Nf = input(‘Nf = ‘); bn = zeros(Nf+1,1); cn = zeros(Nf+1,1); bn(1) = 0; for i = 1:Nf bn(i+1) = (-1)^(i+1)*1/(i*pi); % 计算系数bn cn(i+1) = abs(bn(i+1)); %计算幅度频谱 end t = -5:0.01:5; x = sawtooth(pi*(t+1)); %用sawtooth函数构造周期锯齿脉冲信号 subplot(211); plot(t,x); axis([-5 5 -1.5 1.5]); title(‘周期锯齿脉冲信号‘,‘Fontsize‘,8); xlabel(‘t (单位:s)‘, ‘Fontsize‘,8); subplot(212); k = 0:Nf; stem(k,cn); hold on; plot(k,cn); title(‘幅度频谱‘,‘Fontsize‘,8); xlabel(‘谐波次数‘, ‘Fontsize‘,8); % End
% 周期三角脉冲信号频谱分析 function CTFS_SJ % 绘制周期三角脉冲信号频谱 % Nf：谐波的次数 % an：第1,2,3,...次谐波余弦项系数 display(‘Please input the value of Nf ‘); Nf = input(‘Nf = ‘); an = zeros(Nf+1,1); cn = zeros(Nf+1,1); an(1) = 1/2; an(1) = an(1); for i = 1:Nf an(i+1) = 4*sin(i*pi/2)*sin(i*pi/2)/(i*i*pi*pi); cn(i+1) = abs(an(i+1)); end t = -5:0.01:5; x = (sawtooth(pi*(t+1),0.5)+1)/2; %构造三角脉冲信号 subplot(211); plot(t,x); axis([-5 5 -1.5 1.5]); title(‘周期三角脉冲信号‘,‘Fontsize‘,8); xlabel(‘t (单位:s)‘, ‘Fontsize‘,8); subplot(212); k = 0:Nf; stem(k,cn); hold on; plot(k,cn); title(‘幅度频谱‘,‘Fontsize‘,8); xlabel(‘谐波次数‘, ‘Fontsize‘,8); % End
傅立叶变换：F=fourier(f); F,f应为符号表达式
反傅立叶变换：f=ifourier(F);
代码实现
%ex_4 clear all; R=0.01; t=-2:R:2; f=Heaviside(t+1)-Heaviside(t-1); W1=2*pi*5; N=1000;k=-N:N;W=k*W1/N; F=f*exp(-j*t‘*W)*R; F_r=real(F); figure(1) subplot(2,1,1);plot(t,f); xlabel(‘t‘);ylabel(‘f(t)‘);title(‘f(t)=u(t+1)-u(t-1)‘); subplot(2,1,2);plot(W,F_r); xlabel(‘w‘);ylabel(‘F(w)‘);title(‘f(t)的付氏变换F(w)‘); F_A=abs(F);%幅频特性 F_P=angle(F);%相频特性 figure(2) subplot(2,1,1),plot(W,F_A),xlabel(‘w‘);ylabel(‘abs(F(w))‘);title(‘f(t)幅频特性)‘); subplot(2,1,2),plot(W,F_P),xlabel(‘w‘);ylabel(‘angle(F(w))‘);title(‘f(t)相频特性)‘);
标签：
原文地址：http://www.cnblogs.com/BlueMountain-HaggenDazs/p/4528779.html
