matlab（t）傅⾥叶,MATLAB实现傅⾥叶变换-光⼦的⽇志-⽹易
博客
matlab傅里叶变换的幅度谱和相位谱%% 3.1
%3.1--求图3-60的三⾓波的傅⾥叶级数，并⽤MATLAB，求:
%(1).画出双边幅度谱和相位谱
%(2).若作出 N=3,5,9 时候的波形
%因为3-60的图三⾓波有两个 我选了第⼀个波形 并因为其并未确定 幅度和周期的确切值;
%所以我将这 振幅定义A=4 周期定义 T=4;
clear all
T=4;dt=0.01;
A=4;
t1=-T/2:dt:T/2;
ft1=A/(T/2)*abs(t1);
t=[t1-T t1 t1+T];%⽤于得到三个区间;
ft=repmat(ft1,1,3);%进⾏矩阵复制排列得到与t矩阵匹配的ft
subplot(3,1,1);
plot(t,ft);
xlabel('t');
title('三⾓脉冲的部分图像');
%幅度谱
w0=2*pi/T;
N=3;%可以通过修改这个值 进⾏改变 N 的取值.
K=-N:N;
for k=0:N
Fn(k+1)=quadv([num2str(A),'/2*abs(t)*exp(-i*',num2str(k),'*',num2str(w0),'*t)'],-T/2,T/2)/T;
end
l=length(Fn);%计算出Fn的长度⽤以确定双边谱的长度
F=zeros(1,2*l-1);
F(1:l)=flipud(Fn');F(l:end)=Fn';%对 幅度谱 进⾏双边拓展 成为 F.
subplot(3,2,3);
stem(K*w0,abs(F));% 作出双边幅度谱
%相位谱
ph=angle(F);%求出相位
subplot(3,2,4);
stem(K*w0,ph);
xlabel('nw0');
title('相位谱图像');
%傅⾥叶级数 近似模拟得到三⾓脉冲
K=K';%这⾥的K就是 [-N:N]
ft=F*exp(j*w0*K*t);%从 -N:N 的傅⾥叶级数
subplot(3,1,3);
plot(t,ft);
title('傅⾥叶级数近似表⽰');
%上⾯我做了N=3 的时候的波形 对题设中的 N=5 和N=9 因为只需要修改第20⾏的 N 值就可以得到%% 3.2
%设 f(t)=e^(-2|t|) 求f(t),f(t-2)的傅⾥叶变换，分别作出其幅度谱和相位谱
clear all;
x=15;
T=17;
t0=2;
N=500;
w=linspace(-x,x,N);%对w的取值进⾏线性分割
F=zeros(1,N);%对⽤于存储数据的F进⾏赋初值
for k=1:N
factor=['exp(-j*t*',num2str(w(k)),')'];
F(k)=quad(['exp(-2*abs(t)).*',factor],-T,T);
end
subplot(4,1,1);
plot(w,F);
xlabel('w');
title('f(t)的幅度谱');
subplot(4,1,2);
plot(w,angle(F));
for k=1:N %根据时移特性 f(t-t0)F(jw)e^(-jwt0)进⾏逐项运算得到时移后的波形
F(k)=F(k)*exp(-j*w(k)*t0);
end
subplot(4,1,3);
plot(w,F);
xlabel('w');
title('f(t-2)的幅度谱');
subplot(4,1,4);
plot(w,angle(F));
xlabel('w');
title('f(t-2)的的相位谱');
%% 3.3
%⽤freqs画出下列系统的幅频特性，并确定其是否具有低通.⾼通.或带通特性.
%(1)H(jw)=4/((jw)^3+4(jw)^2+8jw+8) 原题设中的第⼆项为 4(jw)我觉得是少些了个指数所以我就加上了%由输⼊可知 微分⽅程会是: y'''(t)+4y''(t)+8y'(t)+8y(t)=4f(t)
clear all;
a=[1 4 8 8];
b=4;
fs=0.01;
K=4;
w=0:fs:K*pi;
H1=freqs(b,a,w);
subplot(2,1,1);
plot(w,abs(H1));
xlabel('w');
title('H(jw)=4/((jw)^3+4(jw)^2+8jw+8)波形');
%(2)y''(t)+sqrt(2)*y'(t)+y(t)=f''(t)
c=[1 0 0];
d=[1 sqrt(2) 1];
H2=freqs(c,d,w);
subplot(2,1,2);
plot(w,abs(H2));
%信号 f1(t)如图3-90所⽰
%(1) 画出f(t)=f1(t)*cos(50t)的波形.
%(2) 某系统的频率响应为
%    H(jw)=10^4/((jw)^2+26.131(jw)^3+3.4142*10^2(jw)^2+2.6131*10^3(jw)+10^4) %    画出H(jw)的幅频特性和相频特性.
%(3) 将f(t)通过上述系统，画出f(t)和输出信号的幅度谱.
%(4) 画出输出信号的波形.
%(1)
clear all;
dt=0.01;
K=5;
N=500;%修改这⾥可以改变采样点的数⽬
tt=-2:dt:4;
f1=tripuls(tt-1,2);%f1(t)的波形
subplot(4,2,1);
plot(tt,f1);
xlabel('t');
title('f1(t) 的波形');
f=tripuls(tt-1,2).*cos(50*tt);%f(t)的波形
subplot(4,2,2);
plot(tt,f);
xlabel('t');
title('f(t) 的波形');
%(2)由H(jw)易得到输⼊与输出的关系
b=10^4;
a=[1 26.131 3.4142*10^2 2.6131*10^3 10^4];
T=K+1;
w=linspace(-K*pi,K*pi,N);
H=freqs(b,a,w);
Ha=abs(H);%幅频特性
Hph=angle(H);%相频特性
title('Ｈ(jw)的幅频特性');
subplot(4,2,4);
plot(w,Hph);
xlabel('w');
title('H(jw)的相频特性');
%(3)
F=zeros(1,N);
for k=1:N
factor=['exp(-j*t*',num2str(w(k)),')'];
F(k)=quad(['tripuls(t-1,2).*cos(50*t).*',factor],-T,T); Y(k)=F(k)*H(k);
end
subplot(4,2,5);
plot(w,F);
xlabel('w');
title('f(t)的幅度谱');
subplot(4,2,6);
plot(w,Y);
xlabel('w');
title('y(t)的幅度谱');
%(4)
t=linspace(-K*pi,K*pi,N);
for k=1:N  %对输出信号进⾏求算
factor=['exp(j*w*',num2str(t(k)),')'];
y(k)=quad([num2str(Y(k)),'*',factor],-T,T)/(2*pi); end
subplot(4,1,4);
plot(t,y);
xlabel('t');
title('y(t)的波形');
%% 关于page.100 程序的修改

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。