实验四 信号的傅里叶变换

实验内容

题目一

给定一连续 LTI 系统，描述其输入与输出关系的微分方程为：

$$\frac{𝑑 ^2 (𝑦(𝑡))}{𝑑𝑡 ^2} + 5 \frac{𝑑(𝑦(𝑡))}{ 𝑑𝑡} + 3𝑦(𝑡) = 𝑥(𝑡)$$

分别绘制系统的幅频响应、相频响应、频率响应的实部和频率响应的虚部（共四张图），要求频域上计算点数不少于 500 点，各个图标注清晰。其中，求解频率响应可以通过 freqs 函数实现；计算 复数幅度应使用 abs 函数，计算复数相角应使用 angle 函数，计算 复数实部应使用 real 函数，计算复数虚部应使用 imag 函数。

代码

close all;
clear all;
a = [1 5 3];
b = [1];
w = -2 * pi : 0.01 : 2 * pi;
h = freqs(b, a, w);
%幅频响应、相频响应、频率响应的实部和频率响应的虚部（共四张图）
ax1 = subplot(2, 2, 1);
plot(w / pi, abs(h));grid;
title('幅频响应');
xlabel('\omega/\pi');
ax3 = subplot(2, 2, 3);
plot(w / pi, angle(h));grid;
title('相频响应');
xlabel('\omega/\pi');
ax2 = subplot(2, 2, 2);
plot(w / pi, real(h));grid;
title('频率响应的实部');
xlabel('\omega/\pi');
ax4 = subplot(2, 2, 4);
plot(w / pi, imag(h));grid;
title('频率响应的虚部');
xlabel('\omega/\pi');

仿真图

题目二

使用数值以及符号计算方法计算以下信号的傅里叶变换，并分别绘制对应的幅度、相位谱：

$$𝑥_1 (𝑡) = 𝑒 ^{−2|𝑡|} 𝑠𝑖𝑛 4𝑡 \\ 𝑥_2 (𝑡) = 𝑢(𝑡 + 3) − 𝑢(𝑡 − 1)$$

检查要求：两个信号，两种方法，对应的幅度谱和相位谱，一共 8 幅图，请标注清晰，各个信号的四幅图在同一个 figure 中。

代码

close all;
clear 
%符号法
syms t
x1_t = exp(-2*t) *sin(4*t)*heaviside(t)+ exp(2*t) *sin(4*t)*heaviside(-t);
X1 = fourier(x1_t);
x2_t=heaviside(t+3)-heaviside(t-1);
X2=fourier(x2_t);
w=[-10,10];
%数值法
dt = 0.001;
dw = 0.01;
t1 = -10: dt: 10;
w1 = -10:dw:10;
x11_t =  exp(-2*t1).*sin(4*t1).*heaviside(t1)+ exp(2*t1).*sin(4*t1).*heaviside(-t1);
X11 = x11_t * exp(-1j * t1.' * w1) * dt;
x21_t =  heaviside(t1+3)-heaviside(t1-1);
X21 = x21_t * exp(-1j * t1.' * w1) * dt;

%x1图像
figure(1)
subplot(2,2,1);
fplot (abs(X1),w);
title('符号法作幅度谱');
xlabel('\omega');ylabel('|X1(jw)|');grid
ylim([0, 0.5])
subplot(2,2,2);
fplot (angle(X1),w);
title('符号法作相位谱');
xlabel('\omega');ylabel('<X1(jw)');grid
ylim([-2, 2])

subplot(2,2,3);
plot (w1,abs(X11));
title('数值法作幅度谱');
xlabel('\omega');ylabel('|X1(jw)|');grid
ylim([0, 0.5])
subplot(2,2,4);
plot (w1,angle(X11));
title('数值法作相位谱');
xlabel('\omega');ylabel('<X1(jw)');grid
ylim([-2, 2])

%x2图像
figure(2)
subplot(2,2,1);
fplot (abs(X2),w);
title('符号法作幅度谱');
xlabel('\omega');ylabel('|X2(jw)|');grid
subplot(2,2,2);
fplot (angle(X2),w);
title('符号法作相位谱');
xlabel('\omega');ylabel('<X2(jw)');grid

subplot(2,2,3);
plot (w1,abs(X21));
title('数值法作幅度谱');
xlabel('\omega');ylabel('|X2(jw)|');grid
subplot(2,2,4);
plot (w1,angle(X21));
title('数值法作相位谱');
xlabel('\omega');ylabel('<X2(jw)');grid

仿真图

figure1-$x_1(t)$两种方法的图像

figure2-$x_2(t)$两种方法的图像

题目三

验证傅里叶变换的性质：

**线性性质：**计算$5𝑥_1(𝑡) + 𝑥_2(𝑡)$的傅里叶变换（直接计算与使用傅 里叶变换性质计算），分别使用两种方法绘制幅度谱和相位谱，并进行对比；

**卷积性质：**计算$𝑥_1(𝑡) ∗ 𝑥_2(𝑡)$的傅里叶变换（直接计算与使用傅里 叶变换性质计算），分别使用两种方法绘制幅度谱和相位谱，并进 行对比。

傅里叶变换的线性性质体现为：

$$x_1(t)\stackrel{\mathcal{F }}{\leftrightarrow} X_1(jw),x_2(t)\stackrel{\mathcal{F }}{\leftrightarrow} X_2(jw)\\ ax_1(t)+bx_2(t)\stackrel{\mathcal{F }}{\leftrightarrow} aX_1(jw)+bX_2(jw)$$

卷积性质：

$$x_1(t)\stackrel{\mathcal{F }}{\leftrightarrow} X_1(jw),x_2(t)\stackrel{\mathcal{F }}{\leftrightarrow} X_2(jw)\\ x_1(t)∗x_2(t)\stackrel{\mathcal{F }}{\leftrightarrow} X_1(jw)∗X_2(jw)$$

代码

close all;
clear 

dt=0.001;
t = -10: dt: 10;
dw = 0.1;
w = -10:dw:10;
x1_t = exp(-2*t) .*sin(4*t).*heaviside(t)+ exp(2*t) .*sin(4*t).*heaviside(-t);
X1 = x1_t * exp(-1j * t.' * w) * dt;

x2_t=heaviside(t+3)-heaviside(t-1);
X2=x2_t * exp(-1j * t.' * w) * dt;

y_t=5*x1_t+x2_t;
Y=y_t * exp(-1j * t.' * w) * dt;
X3=5*X1+X2;

t1 = -20: dt: 20;
w1=-20:dw:20;
f_t=conv(x1_t,x2_t)*dt;
F=f_t* exp(-1j * t1.' * w1) * dt;
G=X1.*X2;

%线性图像
figure(1)
subplot(2,2,1);
plot (w,abs(X3));
title('5*X_1(jw)+X_2(jw)的幅度谱');
xlabel('\omega');ylabel('|X3(jw)|');grid
ylim([0,4]);
subplot(2,2,2);
plot (w,angle(X3));
title('5*X_1(jw)+X_2(jw)的相位谱');
xlabel('\omega');ylabel('<X3(jw)');grid


subplot(2,2,3);
plot (w,abs(Y));
title('fourier(5x_1(t)+x_2(t))的幅度谱');
xlabel('\omega');ylabel('|Y(jw)|');grid
ylim([0,4]);
subplot(2,2,4);
plot (w,angle(Y));
title('fourier(5x_1(t)+x_2(t))的相位谱');
xlabel('\omega');ylabel('<Y(jw)');grid


%卷积图像
figure(2)
subplot(2,2,1);
plot (w,abs(G));
title('X_1(jw)X_2(jw)的幅度谱');
xlabel('\omega');ylabel('|G(jw)|');grid
subplot(2,2,2);
plot (w,angle(G));
title('X_1(jw)X_2(jw)的相位谱');
xlabel('\omega');ylabel('<G(jw)');grid

subplot(2,2,3);
plot (w1,abs(F));
title('fourier(x_1(t)*x_2(t))的幅度谱');
xlabel('\omega');ylabel('|F(jw)|');xlim([-10,10]);grid
subplot(2,2,4);
plot (w1,angle(F));
title('fourier(x_1(t)*x_2(t))的相位谱');
xlabel('\omega');ylabel('<F(jw)');xlim([-10,10]);grid

仿真图

题目四

使用数值方法计算离散时间傅里叶变换：计算信号$𝑥[𝑛] = ( \frac{1 }{2} ) ^{|𝑛|}$ 的傅里叶变换$ X(𝑒^{𝑗ω})$，绘制幅度谱与相位谱；

$x[n]$是无限长的非周期序列，$X(e^{jw})$是周期的，对其在$[0,\pi]$上求值

$$X(e^{jw})=\sum_{-\infty}^{\infty} x[n]e^{-jwn}=\sum_{0}^{\infty} 0.5^ne^{-jwn}+\sum_{-\infty}^{0} 0.5^{-n}e^{-jwn}\\ =\sum_{n=0}^{\infty} (0.5^ne^{-jwn}+0.5^ne^{jwn})=\frac{1-0.5^2}{1+0.5^2-cosw}$$

代码

clc
clear
close all

dw=0.01; 
w = -2*pi:dw:2*pi; 
X =0.75*1./( 1.25 -cos(w));
magX = abs(X);
angX = angle(X);

subplot(2,1,1)
plot(w/pi,magX);
grid;
title('幅度谱');
xlabel('frequency in pi units');ylabel('Magnitude');
 
subplot(2,1,2)
plot(w/pi,angX);
grid;
title('相位谱')
xlabel('frequency in pi units');ylabel('Radians');
 

仿真图

思考题

使用数值方法计算信号的傅里叶变换时，误差主要有哪些来源？ 请以$𝑥_2 (𝑡)$为例进行说明。

$$f(𝑡) = 𝑢(𝑡 + 3) − 𝑢(𝑡 − 1)$$

$f(t)$ 的傅里叶变换为:

$$F(j \omega)=\int_{-\infty}^{+\infty} f(t) e^{-j \omega t} d t=\lim _{\Delta t \rightarrow 0} \sum_{n=-\infty}^{+\infty} f(n \Delta t) e^{-j \omega n \Delta t} \Delta t$$

若信号为时限信号, 当 $\Delta t$ 足够小, 项数 $N$ 足够多时, 有:

$$\begin{aligned} &F(j \omega)=\Delta t \sum_{n=0}^{N-1} f\left(t_1+n \Delta t\right) e^{-j \omega\left(t_1+n \Delta t\right)} \\ &=\Delta t *\left[f\left(t_1\right), f\left(t_2\right), \ldots, f\left(t_{2 N+1}\right)\right] *\left[e^{-j \omega t_1}, e^{-j \omega t_2}, \ldots, e^{-j \omega t_{2 N+1}}\right]^T \end{aligned}$$

其中 $[\cdot]^T$ 表示矩阵转置。
用 MATLAB 代码表示 $F$ :

$$F=f * \exp \left(-1 j * t .^{\prime} * \omega\right) * d t$$

其中 $f=\left[f\left(t_1\right), f\left(t_2\right), \ldots, f\left(t_{2 N+1}\right)\right]$,$ t=\left[t_1, t_2, \ldots, t_{2 N+1}\right]=[-10,10+dt,10+2*dt,…,10]$,$ \omega=\left[\omega_1, \omega_2, \ldots\right]$。$ dt$ 是 $t$ 的间隔， $t .^{\prime}$ 表示 $t$ 的转置。

数值方法计算傅里叶变换的实质是用足够小的 $\Delta t$和足够多的项数$N$进行求和模拟积分过程，实际过程中我们取定的==$dt$没有足够小==会造成误差，在==区间中取的点没有达到足够多==也造成了误差。对于==$dw$的取样==也会造成误差，绘图过程是集合对应画图，如果$w$的取点不够多，会造成某些取值处点缺失形成误差，如在一开始取$dw=0.1$时第三题的绘图结果，某些点处均产生明显上移：

在计算信号的傅里叶变换方面，数值计算方法和符号计算方法各有什么优缺点？（提示：可从误差、运行时间、实现难度等方面考虑）

• 数值计算方法需要信号是时限信号，因为计算机只能处理有限大小和有限数量的数，取样间隔的大小会对误差造成较大影响，在实现过程中需要考虑较多，实现难度大

• 符号计算方法因为MatLab提供了函数，调用起来比较简单，实现难度小，误差相比数值法手动取样较小。

• 运行时间测试：使用tic,toc分别计算符号法与数值法运行时间：

  dt	dw	符号法计算时间	数值法计算时间
  0.005	0.1	0.101865 秒	0.072984 秒
  0.001	0.1	0.049130 秒	0.180165 秒
  0.001	0.01	0.052693 秒	0.953460 秒
  0.0001	0.1	0.036236 秒	1.181662 秒
  0.0001	0.01	0.063543 秒	16.822756 秒
  0.00001	0.1	0.100490 秒	26.060411 秒

  可以看到，符号法运行时间不受dt影响且速度较快，而数值法随着取样精度的增加运行时间急剧增加。

实验收获与感想

通过本次实验我学习到：

• 利用MATLAB内置函数freqs计算系统幅频、相频响应的方法;
• 熟练掌握了使用 MATLAB 进行傅里叶变换的两种方法——符号法与数值法，仔细比对不同因素对于数值法误差的影响，研究了两种方法的用法、原理与优缺点；
• 利用MATLAB仿真图像验证傅里叶变换的线性性质与卷积性质

此外，随着对MATLAB这一工具的熟练掌握，我还可以自己探究验证傅里叶变换的其他性质，将一些理论学习中的结论可视化，加深理解，做到理论与实践双管齐下。

实际应用过程中常常用到快速傅里叶变换fft()函数，这一部分可以作为拓展内容进一步学习，通过本次实验的学习，收获良多。