实验四--连续时间傅立叶变换.doc

信号与系统实验报告

实验四 连续时间傅立叶变换
§4.1连续时间傅立叶变换的数字近似
1． 求CTFT的解析表达式。可将看作，。
g=sym('exp(-2*t)*Heaviside(t)');
g2=subs(g,'-t','t');
x=g+g2;
fx=fourier(x);

2． 创建一个向量，它包含了在区间t=[0:tau:T-tau] 上（其中和），信号的样本。
clc;
t=0:0.01:10-0.01;
y=exp(-2*(t-5)).*(0.5+0.5*sign(t-5))+exp(2*(t-5)).*(0.5+0.5*sign(-t+5));
plot(t,y);

3． 键入y=fftshift(tau*fft(y))计算样本。因为对于基本上为零，就能近似用个样本分析中计算出信号的CTFT。
clc;
t=0:0.01:10-0.01;
y=exp(-2*(t-5)).*(0.5+0.5*sign(t-5))+exp(2*(t-5)).*(0.5+0.5*sign(-t+5));
y=fftshift(0.01.*fft(y));
y=abs(y);
plot(t,y);
axis([4,6,-0.1,1.2]);

4.构造一个频率样本向量w，它按照
>> w=-(pi/tau)+(0:N-1)*(2*pi/(N*tau));
与存在向量Y中的值相对应。
5.因为是通过时移与相联系的，所以CTFT就以线性相移项与相联系。利用频率向量w直接由Y计算的样本，并将结果存入x中。
clc;
t=0:0.01:10-0.01;
tau=0.01;
N=10/0.01;
y=exp(-2*(t-5)).*(0.5+0.5*sign(t-5))+exp(2*(t-5)).*(0.5+0.5*sign(-t+5));
y=fftshift(0.01.*fft(y));
w=-(pi/tau)+(0:N-1)*(2*pi/(N*tau));
x=exp(j*5*w).*y;

6.利用abs和angle画出在w标定的频率范围内X的幅值和相位。对于相同的值，也画出在1中所导出的解析式表达式的幅值和相位。CTFT的近似值与解析导得的相符吗?若想在一张对数坐标上画出幅值，可以用semilogy，这是会注意到，在较高的频率上近似不如在较低的频率上好。因为用了样本近似，所以在时间段长度内，信号变化不大的那些信号的频率分量近似程度会更好一些。
clc;
tau=0.01;
T=10;
t=0:tau:T-tau;
N=T/tau;
y=exp(-2*(t-5)).*(0.5+0.5*sign(t-5))+exp(2*(t-5)).*(0.5+0.5*sign(-t+5));
y=fftshift(0.01.*fft(y));
w=-(pi/tau)+(0:N-1)*(2*pi/(N*tau));
x=exp(i*5*w).*y;
xp=abs(x);
xf=angle(x);
subplot(211)
plot(t,xp);
subplot(212)
plot(t,xf);

7.利用abs和angle画出Y的幅值和相位，它们与X的图比较后怎样？能估计到这一结果吗？
clc
syms s w;
g=sym('exp(-2*s)*Heaviside(s)');
g2=subs(g,'-s','s');
y=g+g2;
fw=fourier(y,s,w);
ff=atan(imag(fw)/real(fw));
fp=abs(fw);
tau=0.01;
T=10;
t=0:tau:T-tau;
N=T/tau;
y=exp(-2*(t-5)).*(0.5+0.5*sign(t-5))+exp(2*(t-5)).*(0.5+0.5*sign(-t+5));
y=fftshift(0.01.*fft(y));
w=-(pi/tau)+(0:N-1)*(2*pi/(N*tau));
x=exp(i*5*w).*y;
xp=abs(x);
xf=angle(x);
subplot(211)
hold on
plot(t,xp);
ezplot(fp,-10:10)
hold off
subplot(212)
hold on
plot(t,xf);
ezplot(ff,-10:10)
hold off

§4.2连续时间傅立叶变换性质
1．键入Y=fftshift(fft(y))，计算向量Y的傅立叶变换。键入
>> w=[-pi:2*pi/N:pi-pi/N]*fs;
将对应的频率值存入向量w中。利用w和Y在区间内画出该连续时间傅立叶变换的幅值。
函数ifft是fft的逆运算。对于偶数长度的向量，fftshift就是它本身的逆。对于向量Y，N=8192，这个逆傅立叶变换能用键入以下命令而求得
>> y=ifft(fftshift(Y));
>> y=real(y);
由于原时域信号已知是实的，所以这里用了函数real。然而，在fft和ifft中的数值舍入误差都会在y中引入一个很小的非零虚部分量。一般说来，逆CTFT不必是一个实信号，而虚部可以包含有显著的能量。当已知所得信号一定是实信号时，并且已经证实所除掉的虚部分量是没有意义的，real函数才能用于ifft的输出上
load splat
y=y(1:8192);
N=8192;
fs=8192;
Y=fftshift(fft(y));
sound(y,fs);
w=(-pi:2*pi/N:pi-pi/N)*fs;
subplot(211)
plot(w,Y);
title('Y');
y=ifft(fftshift(Y));
y=real(y);
subplot(212)
plot(w,y);
title('y');

2． 置Y1=conj(Y)并将Y1的逆傅立叶变换存入Y1中，用real(y1)以确保y1是实的，用sound(y1,fs)将y1放出。已知的逆傅立叶变换是如何与联系的，能解释刚才听到的是什么吗？
load splat
y=y(1:8192);
N=8192;
fs=8192;
Y=fftshift(fft(y));
Y1=conj(Y);
y1=ifft(fftshift(Y1));
sound(y1,fs);
w=(-pi:2*pi/N:pi-pi/N)*fs;
plot(w,y1);

答：刚才听到的是y信号反过来放的声音。

的CTFT可以用它的幅值和相位写成

式中。对于许多信号，单独用相位或幅值都能构造出一个有用的信号的近似。例如，考虑信号和，其CTFT为

3.只要是实信号，用解析方法说明和一定是实的。
解：因为即y(t)=※又因为是实信号故和一定是实的。
4.构造一个向量Y2等于Y的幅值，并将Y2的逆傅立叶变换存入向量y2中，用sound放出这个向量。
load splat
y=y(1:8192);
N=8192;
fs=8192;
Y=fftshift(fft(y));
sound(y,fs);
w=(-pi:2*pi/N:pi-pi/N)*fs;
Y2=abs(Y);
y2=ifft(fftshift(Y2));
sound(y2,fs);
plot(w,y2);

5.构造一个向量Y3，它有与Y相同的相位，但是幅值对每个频率都等于1。并将Y3的逆傅立叶变换存入向量y3中，用sound放出这个向量。
load splat
y=y(1:8192);
N=8192;
fs=8192;
Y=fftshift(fft(y));
sound(y,fs);
w=(-pi:2*pi/N:pi-pi/N)*fs;
Y3=Y./abs(Y);
y3=ifft(fftshift(Y3));
sound(y3,fs);
plot(w,y3);

6.根据刚才听到的这两个信号，代表一个声音信号你认为傅立叶变换的那个部分是最关键的：幅值或相位？
答：相位是最关键的。
7.用向量y创建一个向量y4，它包含有本该以8192Hz从采样所得到的样本。
load splat
y=y(1:8192);
N=8192;
fs=8192;
sound(y,fs);
y4=y(1:2:8192);
sound(y4,fs);
w=(-pi:2*pi/N:pi-pi/N)*fs;
w4=(-pi:2*pi/(N/2):pi-pi/(N/2))*(fs/2);
subplot(211)
plot(w,y);
title('y');
subplot(212)
plot(w4,y4);
title('y4');

8.用y4=sound(y4,fs)放出y4。利用比较y4的傅立叶变换与y的傅立叶变换，能说明在高音上的变化吗？信号压缩是如何影响它的傅立叶变换的？
load splat
y=y(1:8192);
N=8192;
fs=8192;
sound(y,fs);
y4=y(1:2:8192);
sound(y4,fs);
Y=fftshift(fft(y));
Y4=fftshift(fft(y4));
w=(-pi:2*pi/N:pi-pi/N)*fs;
w4=(-pi:2*pi/(N/2):pi-pi/(N/2))*(fs/2);
subplot(211)
plot(w,Y);
title('Y');
subplot(212)
plot(w4,Y4);
title('Y4');

信号在时域中压缩（a>1）等效于在频域中扩展。
9．创建向量x，，它由下式给出

注意，x是一个长度为2*N的向量。
load splat
y=y(1:8192);
N=8192;
fs=8192;
sound(y,fs);
j=1;
for i=1:2:(2*8192)
x(i)=0;
x(i+1)=y(j);
j=j+1;
end
10.利用函数filter完成在x上的线性内插。这里要用到的线性内插器的单位冲激响应是h=[1 2 1]/2。
load splat
y=y(1:8192);
N=8192;
fs=8192;
sound(y,fs);
j=1;
for i=1:2:(2*8192)
x(i)=0;
x(i+1)=y(j);
j=j+1;
end
h=[1 2 1];
a=2;
y=filter(h,a,x);
11.用sound(y5,fs)放出y5。用比较y5和y的傅立叶变换，能解释在音调上的变化吗？
load splat
y=y(1:8192);
N=8192;
fs=8192;
sound(y,fs);
j=1;
w=(-pi:2*pi/N:pi-pi/N)*fs;
w5=(-pi:2*pi/(N*2):pi-pi/(N*2))*(fs*2);
for i=1:2:(2*8192)
y5(i)=0;
y5(i+1)=y(j);
j=j+1;
end
sound(x,fs);
Y=fftshift(fft(y));
Y5=fftshift(fft(y5));
subplot(211)
plot(w,Y);
subplot(212)
plot(w5,Y5)

§4.3连续时间傅立叶变换的符号计算
1．定义符号表达式x1和x2代表下面连续时间信号：

需要用函数Heaviside来表示单位阶跃函数。
x1=sym('Heaviside(t)*0.5*exp(-2*t)');
x2=sym('Heaviside(t)*exp(-4*t)');
2． 对于1中所定义的和，用解析方法计算它们的CTFT在的值，即（不应该先求来作这道题）CTFT在的值是怎样与时域信号关联的？
clc;
syms t
x1=sym('Heaviside(t)*0.5*exp(-2*t)');
x2=sym('Heaviside(t)*exp(-4*t)');
fx1=simple(int(x1, t, -inf, inf));
fx2=simple(int(x2, t, -inf, inf));
3．1所定义的信号中，哪一个在时域衰减得更快？根据这一点，你能预期在频域哪一个衰减得更快？
x1=sym('Heaviside(t)*0.5*exp(-2*t)');
x2=sym('Heaviside(t)*exp(-4*t)');
subplot(211);
ezplot(x1,-1:5);
axis([-1,5,-0.1,0.8])
subplot(212);
ezplot(x2,-1:5);
axis([-1,5,-0.1,1.2])

X2在时域衰减得更快,X1在频域衰减得更快.

4． 用函数fourier计算和得CTFT。定义x1和x2是由fourier产生的符号表达式。用ezplot产生和的幅值图。这些图能对2和3中的答案进行确认吗？
syms t w;
x1=sym('Heaviside(t)*0.5*exp(-2*t)');
x2=sym('Heaviside(t)*exp(-4*t)');
fw1=fourier(x1,t,w);
fw2=fourier(x2,t,w);
ffw1=maple('convert',fw1,'piecewise');
ffw2=maple('convert',fw2,'piecewise');
ffp1=abs(ffw1);
ffp2=abs(ffw2);
subplot(211);
ezplot(ffp1,-100:100);
axis([-20,20,-0.1,0.3]);
title('fourier(x1)');
subplot(212);
ezplot(ffp2,-100:100);
axis([-20,20,-0.1,0.3]);
title('fourier(x2)');

5．定义符号表达式y1代表下面连续时间信号：

它可以作为两个Heaviside函数之差。
y1=sym('Heaviside(t+2)-Heaviside(t-2)');
6.用解析方法求的CTFT，。
syms t w;
y=sym('Heaviside(t+2)-Heaviside(t-2)');
F=int(y*exp(-1i*w*t),t, -inf, inf);
ezplot(F,-4*pi:4*pi);
axis([-4*pi 4*pi 0 5])
title('fourier(y1)');

7．定义符号表达式y2表示信号。你能像对y1那样用两个Heaviside函数之差来完成，或者恰当地对y1应用subs。
syms t w;
y1=sym('Heaviside(t+2)-Heaviside(t-2)');
y2=subs(y1,'t-2','t');
fw=fourier(y2,t,w);
ffw=maple('convert',fw,'piecewise');
ffp=abs(ffw);
ezplot(ffp,-4*pi:4*pi);
axis([-4*pi 4*pi 0 2.5])
title('fourier(y2)');

8．利用fourier求y1和y2的CTFT，并将它们存入Y1和Y2中。倘若Y1不是你所期望得到的表达式，那么试试在所得表达式上用simple以便得出更为熟悉的形式。
syms t w;
y1=sym('Heaviside(t+2)-Heaviside(t-2)');
y2=subs(y1,'t-2','t');
Y1=fourier(y1,t,w);
Y2=fourier(y2,t,w);

9．用ezplot产生和的幅值图。比较这两张图情况如何？由这两个信号在时域之间的关系能预测到这个结果吗？
syms t w;

y1=sym('Heaviside(t+2)-Heaviside(t-2)');
y2=subs(y1,'t-2','t');
Y1=fourier(y1,t,w);
fY1=maple('convert',Y1,'piecewise');
fpY1=abs(fY1);

Y2=fourier(y2,t,w);
fY2=maple('convert',Y2,'piecewise');
fpY2=abs(fY2);

hold on
ezplot(fpY1,-4*pi:4*pi);
ezplot(fpY2,-4*pi:4*pi);
axis([-4*pi 4*pi 0 4.5])

10． 下面几部分的CTFT。将写成和两个信号之和。将选为因果信号，选为反因果信号，即。用解析方法计算的CTFT，。
syms t w;
v1=sym('exp(-2*t)*Heaviside(t)');
v2=sym('exp(2*t)*Heaviside(-t)');
v=v1+v2;
vf=int(v*exp(-i*w*t),t, -inf, inf);
ezplot(vf,-10:10);

11． 用fourier求v的CTFT的符号表达式V。这个表达式等效于在10中用解析法求得的表达式吗？
syms t w;
v1=sym('exp(-2*t)*Heaviside(t)');
v2=sym('exp(2*t)*Heaviside(-t)');
v=v1+v2;
vf=fourier(v,t,w);
vw=maple('convert',vf,'piecewise');
vp=abs(vw);
ezplot(vp,-10:10);

12． 定义f是信号的符号表达式，用fourier定义F是f的CTFT的符号表达式。注意，F中含有一个未被求值的积分。这个未被求值的积分对所有值都收敛吗？
syms t a w;
f=sym('exp(-a*t)*Heaviside(t)');
F=fourier(f,t,w);
当a=0时此未被求值的积分不收敛。
13．用subs设置F1中的值等于5，然后在这个置换的结果应用simple。所得结果是所期望的吗？
syms t a w;
f=sym('exp(-a*t)*Heaviside(t)');
F=fourier(f,t,w);
F1=simple(subs(F,5,'a'));