周期三角波傅里叶级数例题_MATLAB信号与系统分析(五)——连续时间信号的频谱分析...

一、实验目的:

1、掌握傅立叶级数(FS),学会分析连续时间周期信号的频谱分析及MATLAB实现;

2、掌握傅立叶变换(FT),了解傅立叶变换的性质以及MATLAB实现。

二、利用符号运算求傅里叶级数的系数

1、复习几个函数:

F1=int(f,v,a,b) — 对f表达式的v变量在(a,b)区间求定积分

F2=subs(s,OLD,NEW)-用新变量NEW代替S中的指定变量OLD。

F3=vpa(x,n) : 显示可变精度计算;x为符号变量,n表示要精确计算的位数。

2、周期函数的傅里叶级数的形式

3、利用符号运算求傅立叶级数的系数

代码抄袭如下:

%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)/2Nf=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.01n= 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

4、周期信号的频谱分析:

%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)/2A(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('连续时间函数的双边谱')

5、利用MATLAB实现典型周期信号的频谱

(1)周期方波脉冲频谱的Matlab实现

%周期方波信号频谱分析

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

(2)周期锯齿波脉冲频谱的Matlab实现

%周期锯齿脉冲信号频谱分析

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

(3)周期三角波脉冲频谱的Matlab实现

%周期三角脉冲信号频谱分析

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

三、非周期函数的傅立叶变换

1、利用符号函数求傅立叶变换

傅立叶变换:F=fourier(f); F,f应为符号表达式

反傅立叶变换:f=ifourier(F);

2、连续时间信号傅立叶变换的数值计算

代码实现

%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)相频特性)');

相关资源:MATLAB傅里叶级数
表情包
插入表情
评论将由博主筛选后显示,对所有人可见 | 还能输入1000个字符
相关推荐
©️2020 CSDN 皮肤主题: 1024 设计师:白松林 返回首页