ctfs.m

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

%========================================================================
%               周期信号的傅里叶级数   程序名：CTFS
%   本程序根据三角形式傅里叶级数的系数的计算公式，采用符号计算法计算出信号
%   的傅里叶级数的系数，然后再根据你输入的谐波次数合成（逼近）原周期信号。
% 	Nf  	合成信号的谐波的次数,由实验者从键盘输入，谐波次数可以不限，但
%   建议取50以内，因为次数太大，计算机的运算量越大，所需计算时间越长。
%	A_sym	第1元素是直流项，其后元素依次是1,2,3...次谐波cos项展开系数
%	B_sym	第2,3,4,...元素依次是1,2,3...次谐波sin项展开系数
%   tao     脉冲宽度。周期矩形信号的宽度为0.5s，锯齿波为4s，三角波为2s。
%   T=4     信号周期为4s
%========================================================================
clear, close all
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('2*(t)/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
Fn=int(x*exp(-j*2*pi*n*t/T)/T,t,-T/2,T/2);
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码数值数组
    F_sym(k+1)=double(vpa(subs(Fn,n,k)));
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项展开系数
Fn=F_sym
Fn(1)=a_n(1);%A0%(1);
%===========================================================================
%                 下段程序计算原周期信号的时域波形
%===========================================================================
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.*(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('Approximation of aperiodic signal')
%axis([-8,8,-1.2,1.2]);
%====================================
%      下段程序绘制信号的频谱图
%====================================
figure(2);
n=0:Nf;
c_n=sqrt(a_n.^2+b_n.^2);
subplot(2,1,1);
stem(n,abs(Fn),'.r');
%stem(n,c_n,'.r');
title('Amplitude spectra of aperiodic signal |cn|');
%title('Magnitude spectra of aperiodic signal  Fn');
fi=-atan(abs(imag(Fn))./(abs(real(Fn))+eps));
for L=1:length(n);
    
if b_n(L)==0
    if a_n(L)<=-0.0001
        fi(L)=pi;
    else
        fi(L)=0;
    end
end
end
subplot(2,1,2);
stem(n,fi,'.k');
title('Phase spectra of aperiodic signal');
xlabel('index n');
%====================================