运动裤裤裆老是磨坏:MATLAB 实现傅里叶变换

MATLAB 实现傅里叶变换

信号处理和分析、傅里叶变换 2010-03-25 21:32:27 阅读113 评论1   字号：大中小 订阅

%% 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));% 作出双边幅度谱

xlabel('nw0');

title('双边幅度频谱');

%相位谱

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));

xlabel('w');

title('f(t)的相位谱');

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));

xlabel('w');

title('y``(t)+sqrt(2)*y`(t)+y(t)=f``(t)波形');

%% 3.4

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

subplot(4,2,3);

plot(w,Ha);

xlabel('w');

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 程序的修改

%以下是我修改的程序

clear all;

T=4;

width=2;

A=0.5;

t1=-T/2:0.01:T/2;

ft1=0.5*[abs(t1)

t=[t1-T t1 t1+T];

ft=repmat(ft1,1,3);

subplot(3,1,1);

plot(t,ft);

xlabel('t');

title('original square waveform');

w0=2*pi/T;

N=6;

K=0:N;

for k=0:N

    Fn(k+1)=quadv([num2str(A),'*rectpuls(t,2)*exp(-j*',num2str(k),'*',num2str(w0),'*t)'],-T/2,T/2)/T;

   %Fn(k+1)=quadv([num2str(A),'*rectpuls(t,2)*exp(-j*',num2str(k),'*'num2str(w0)'*t')],-T/2,T/2)/T;

   %注释中的是书上的程序,可是书写上有问题.在num2str(w0)这里逗号少写了两个

end

subplot(3,2,3);

stem(K*w0,abs(2*Fn));%这里的单边谱的Fn应该要乘以2

xlabel('nw0');

title('magnitude spectrum');

%phase

ph=angle(Fn);

subplot(3,2,4);

stem(K*w0,ph);

xlabel('nw0');

title('phase spectrum');

%sythesis

L=length(Fn);

F=zeros(1,2*L-1);

F(1:L)=flipud(Fn');F(L:end)=Fn';

K=[-N:N]';%f(t)是应该从-N:N

ft=F*exp(j*w0*K*t);

subplot(3,1,3);

plot(t,ft);

title('sythesized square wavform');