关于MATLAB实现离散时间傅里叶变换（DTFT）——傅里叶变换的学习与认识（一）

傅里叶变换

法国工程师傅里叶指出，一个“任意”的周期函数x(t)可以分解为无穷多个不同频率的正弦信号的和，这既是傅里叶级数。求解傅里叶傅里叶系数的过程就是傅里叶变换。傅里叶级数（针对周期信号，要求信号在一个周期内的能量是有限的）和傅里叶变换（针对非周期信号，要求信号在整个时间区内的能量是有限的）又统称为傅里叶分析或谐波分析

1.连续时间的傅里叶变换

其目的是以正弦函数(正弦和余弦函数统称为正弦函数)或虚指数
$e^{jwt}$
为基本信号，将任意的连续时间信号表示为一系列不同频率的正弦函数或虚指数函数之和（对于周期信号）或积分（对于非周期信号）

2.离散时间信号的傅里叶变换（DTFT）

它是以虚指数函数
$e^{j(\frac{2\pi}{N})kn}$
或
$e^{j\theta k}$
为基本信号，将任意离散时间信号表示为N个不同频率指数的虚指数之和（对于周期信号）或积分（对于非周期信号）。

同时，为了实现用数字计算机进行傅里叶变换的计算，定义了离散傅里叶变换（DFT）

序列x[n]的离散时间傅里叶变换X(e^jw)是w的连续函数。由于数据在 MATLAB中以向量的形式存在,X(e^jw)只能在一个给定的离散频率的集合中计算。只有类似
$X(e^{jw})= \frac{p_0+p_1 e^{-jw}+...+p_M e^{-jwM}} {d_0+d_1 e^{-jw}+...+d_N e^{-jwN}}$
形式的e^(-jw)的有理函数,才能计算其离散时间傅里叶变换。

在下面的练习中,我们将学习怎样用 MATLAB来求取并画出序列的离散时间傅里叶变换

例.计算离散时间傅里叶变换

如前面提到的式子所示,序列x[n]的离散时间傅里叶变换X(e^(jw)),可以用 MATLAB函数 freqz非常方便地在给定的L个离散频率点w=w_L处进行计算。由于X(e^jw)是w的连续函数,需要尽可能大地选取L的值,使得产生的图形和真实离散傅里叶变换的图形尽可能一致。在 MATLAB中, fretz计算出序列{p0 p1… pM}和序列{d0 d1…dN}的L点离散傅里叶变换,然后对其离散傅里叶变换值相比得到X(e^jw L),L=1,2,…,L。为更加方便快速地运算,应将L的值选为2的幂,如256或者512。

例：离散时间傅里叶变换的原始序列为
$H(e^{jw})=\frac{2+e^{-jw}}{1-0.6e^{-jw}}$
实现代码

% Program P3_1
% 求取信号的离散时间傅里叶变换
clf;
% 计算离散时间傅里叶变换的频率样本
w = -4*pi:8*pi/511:4*pi;  %频率点
num = [2 1];  %系统传递函数的分子的系数
den = [1 -0.6];  %系统传递函数的分母的系数
h = freqz(num, den, w);%h为求取的滤波器频率响应函数
% 画出DTFT
subplot(2,2,1)
plot(w/pi,real(h));grid
title('H(e^{j\omega})的实部')
xlabel('\omega /\pi');
ylabel('振幅');
subplot(2,2,2)
plot(w/pi,imag(h));grid
title(' H(e^{j\omega})的虚部')
xlabel('\omega /\pi');
ylabel('振幅');
pause  %暂停，按任意键继续执行下面的操作
subplot(2,2,3)
plot(w/pi,abs(h));grid
title(' |H(e^{j\omega})|幅度谱')
xlabel('\omega /\pi');
ylabel('振幅');
subplot(2,2,4)
plot(w/pi,angle(h));grid
title('相位谱arg[H(e^{j\omega})]')
xlabel('\omega /\pi');
ylabel('以弧度为单位的相位');

运行结果

例：计算有限长序列的离散时间傅里叶变换

g[n]={1 3 5 7 9 11 13 15 17}

实现代码

% Program P3_1
% 求取有限长序列的离散时间傅里叶变换
clf;
% 计算离散时间傅里叶变换的频率样本
w = -4*pi:8*pi/511:4*pi;  %频率点
a  = [1 3 5 7 9 11 13 15 17 ];  %系统传递函数的分子的系数  
 h = freqz(a, 1, w);%h为求取的滤波器频率响应函数，1为系统传递函数的分母的系数
% 画出DTFT
subplot(2,2,1)
plot(w/pi,real(h));
grid;
title('H(e^{j\omega})的实部')
xlabel('\omega /\pi');
 ylabel('振幅');   
subplot(2,2,2)
plot(w/pi,imag(h));
grid;
title(' H(e^{j\omega})的虚部')
xlabel('\omega /\pi');
ylabel('振幅');
subplot(2,2,3)
plot(w/pi,abs(h));
grid;
title(' |H(e^{j\omega})|幅度谱')
xlabel('\omega /\pi');
ylabel('振幅');
subplot(2,2,4)
plot(w/pi,angle(h));
axis([0 1 -4 4]);
grid;
title('相位谱arg[H(e^{j\omega})]')
xlabel('\omega /\pi');
ylabel('以弧度为单位的相位');

运行结果

对相位谱截取[0,1]之间的图像发现明显的相位谱跳变。相位到一定负值就会跳变为正值，正弦信号的周期为2pi，即-pi到pi。所以，当低于-pi时就会自动加2kpi，发生相位跳变。unwrap命令可移除跳变。

要计算一个系统相频特性，就要用到反正切函数，计算机中反正切函数规定，在一、二象限中的角度为0～pi，三四象限的角度为0～-pi。
若一个角度从0变到2pi，但实际得到的结果是0～pi，再由-pi～0，在w=pi处发生跳变，跳变幅度为2pi，这就叫相位的卷绕。
unwrap( )就是解卷绕，使相位在pi处不发生跳变，从而反应出真实的相位变化 "
查看 unwrap 的帮助文档可以发现 unwrap 还可以输入一个参数 tol，默认tol = pi。也可以根据情况修改 tol。
参考链接：https://blog.csdn.net/zb1165048017/article/details/50381021

实现代码

% 求取有限长序列的离散时间傅里叶变换，
%并对相位谱进行解卷绕
clf;
% 计算离散时间傅里叶变换的频率样本
w = -4*pi:8*pi/511:4*pi;  %频率点
a  = [1 3 5 7 9 11 13 15 17 ];  %系统传递函数的分子的系数  
 h = freqz(a, 1, w);%h为求取的滤波器频率响应函数，1表示系统传递函数的分母的系数为1
% 画出DTFT
subplot(2,2,1)
plot(w/pi,real(h));
grid;
title('H(e^{j\omega})的实部')
xlabel('\omega /\pi');
 ylabel('振幅');   
subplot(2,2,2)
plot(w/pi,imag(h));
grid;
title(' H(e^{j\omega})的虚部')
xlabel('\omega /\pi');
ylabel('振幅');
subplot(2,2,3)
plot(w/pi,abs(h));
grid;
title(' |H(e^{j\omega})|幅度谱')
xlabel('\omega /\pi');
ylabel('振幅');
subplot(2,2,4)
y=unwrap(angle(h));
plot(w/pi,y);
grid;
title('相位谱arg[H(e^{j\omega})]')
xlabel('\omega /\pi');
ylabel('以弧度为单位的相位');

运行结果

未解卷绕之前图像

转载：https://blog.csdn.net/qq_40644291/article/details/102490750