ctfs2.m

来自「经典《信号与系统》教程的matlab例程,对深入理解信号与系统相关概念有很大帮助」· M 代码 · 共 76 行

clear
syms t  n   k  x 
T=4;
Nf=input('请输入信号合成的谐波次数：');
%选择信号类型，确定出信号一个周期内的符号表达式
q=input('请选择信号代码（1-周期矩形波信号，2-周期三角波信号，3-周期锯齿波信号）：');
p=(q~=1)&(q~=2)&(q~=3);
while p==1;
      q=input('请选择信号代码（1-周期矩形波信号，2-周期三角波信号，3-周期锯齿波信号）：');
      p=(q~=1)&(q~=2)&(q~=3);
end
if q==1
    x=sym('Heaviside(t+1)-Heaviside(t-1)');
end
if q==2
    x=sym('(t+2)/4*(Heaviside(t+2)-Heaviside(t-2))');
end
if q==3
    x=sym('(-abs(t)+1)*(Heaviside(t+1)-Heaviside(t-1))');
end
% 下面的程序段计算傅里叶级数的系数
A0=int(x/T,t,-T/2,T/2);                   %求出三角函数展开系数A0
As=int(2*x*cos(2*pi*n*t/T)/T,t,-T/2,T/2); %求出三角函数展开系数As
Bs=int(2*x*sin(2*pi*n*t/T)/T,t,-T/2,T/2); %求出三角函数展开系数Bs
A_sym(1)=double(vpa(A0));                %获取串数组A0所对应的ASC2码数值数组
for k=1:Nf
A_sym(k+1)=double(vpa(subs(As,n,k))); %获取串数组A所对应的ASC2码数值数组
    B_sym(k+1)=double(vpa(subs(Bs,n,k))); %获取串数组B所对应的ASC2码数值数组
end
disp('显示an列表：第1元素是直流项，其后元素依次是1,2,3...次谐波cos项展开系数');
a_n=A_sym %输出cn:第1元素是直流项，其后元素依次是1,2,3...次谐波cos项展开系数
disp('显示bn列表：第2,3,4,...元素依次是1,2,3...次谐波sin项展开系数');
b_n=B_sym  %输出dn:第2,3,4,...元素依次是1,2,3...次谐波sin项展开系数
% 下段程序计算原周期信号的时域波形
t=-8:0.01:8;
f=a_n(1);
for l=1:Nf;
    f=f+a_n(l+1).*cos(2*l*pi*t/T)+b_n(l+1).*sin(2*l*pi*t/T); 
end;
if q==1
    y=0;
   for r=-2:2;
      y=y+u(t+r*4+1)-u(t+r*4-1);
   end;
end
if q==2
    y=0;
   for r=-2:2;
      y=y+(t+r*4+2)/4.*(u(t+r*4+2)-u(t+r*4-2));
   end;
end
if q==3
    y=0;
   for r=-2:2;
      y=y+(-abs(t+r*4)+1).*(u(t+r*4+1)-u(t+r*4-1));
   end;
end
% 下段程序绘制原周期信号和逼近信号的波形图
figure(1);
clf
plot(t,y,'r:');hold on;
plot(t,f);
title('非正弦周期信号形波的合成')
axis([-8,8,-0.4,1.3]);
% 下段程序绘制信号的频谱图
figure(2);
n=0:Nf;
c_n=sqrt(a_n.^2+b_n.^2);
subplot(2,1,1);
stem(n,c_n,'.r');
title('非正弦周期信号的幅度谱|cn|');
fi=-atan(b_n./a_n);
subplot(2,1,2);
stem(n,fi,'.k');
title('非正弦周期信号的相位谱');