Matlab模拟傅里叶变换

傅里叶变换是我们最早开始接触的时频域变换方法,虽然经常使用,知道怎么用纸笔计算,但是还从来没有在电脑中模拟过,正好现在开始学习数字信号处理,借着这个机会再学习如何在电脑上模拟傅里叶变换。

以下大部分内容来自Digital Signal Processing Using Matlab数字信号处理教程 程佩青

此次选择的软件平台为Matlab。

由于Matlab无法处理无限长序列,所以需要处理的信号必须是有限长的。

连续时间傅里叶变换

傅里叶变换的公式为:

为了在计算机中模拟傅里叶变换,我们将积分变为求和的方式,上下限也从正无穷到负无穷变为一段长度M,dt需要尽可能小

在Matlab中,函数的自变量因变量的集合都是使用矩阵来存储的,从矩阵的角度来看傅里叶变换的公式如下:

角频率向量定义为

时间向量定义为

因此矩阵指数可写为

整个傅里叶变换可写为

Xa = xa * exp(-1j*t'*W) * Dt;

具体实现

其实下面这个例子是Digital Signal Processing Using Matlab中的,来自P64页,不过想到都看到这里了还要读者翻书不太好,就一起放上来了。

定义

先进行数学上的分析,

MATLAB实现如下:

% Analog Signal
Dt = 0.00005;
t = -0.005:Dt:0.005;
xa = exp(-1000*abs(t));

% Continuous-time Fourier Transform
Wmax = 2*pi*2000;
K = 500;
k = 0:1:K;
W = k*Wmax/K;

Xa = xa * exp(-1j*t'*W) * Dt;
Xa = abs(Xa);

W = [-fliplr(W), W(2:501)];
Xa = [fliplr(Xa), Xa(2:501)];

subplot(2,1,1); 
plot(t*1000,xa);
xlabel('t in msec.'); 
ylabel('xa(t)');
title('Analog Signal');

subplot(2,1,2); 
plot(W/(2*pi*1000),Xa*1000);
xlabel('Frequency in KHz'); ylabel('Xa(jW)*1000');
title('Continuous-time Fourier Transform');

运行效果如下:

如果想确认变换的正确性,可以在运行完上面这个脚本后,在命令行输入

plot(W/(2*pi*1000),(0.002./(1+(W./1000).^2))*1000);
xlabel('Frequency in KHz'); ylabel('Xa(jW)*1000');

运行效果如下:

这时会发现,根据上面推导的变换公式直接plot出的图形和变换后得到的图形是一样的,这样可以确定变换的正确性。

存在问题

目前存在的问题是,对于复函数的变换结果不正确。我想了很多天都找不出问题所在,只能暂时放弃,等以后有机会再研究。

离散时间傅里叶变换

下面是对上一个例子中的模拟输入信号做离散化,然后再进行离散傅里叶变换。

为了体现Nyquist定理,将使用两种不同的采样频率

  1. 使用Fs=5000sam/sec采样来获得x1(n)

  2. 使用Fs=1000sam/sec采样来获得x2(n)


% Analog Signal
Dt = 0.00005;
t = -0.005:Dt:0.005;
xa = exp(-1000*abs(t));

% Discrete-time Signal 
Ts = 0.0002;
n = -25:1:25;
x = exp(-1000*abs(n*Ts));

% Discrete-time Fourier transform
K = 500;
k = 0:1:K;
w = pi*k/K;

X = x*exp(-j*n'*w); X = real(X);

w = [-fliplr(w), w(2:K+1)];
X = [fliplr(X), X(2:K+1)];

subplot(2,1,1);plot(t*1000,xa);
xlabel('t in msec.'); 
ylabel('x1(n)');
title('Discrete Signal');hold on;

stem(n*Ts*1000,real(x));gtext('Ts=0.2 msec');hold off;

subplot(2,1,2);plot(w/pi,X);
xlabel('Frequency in pi units');ylabel('X1(w)');
title('Discrete-time Fourier Transform');

Fs=5000sam/sec

xa(t)的频率为2KHz,因此它的Nyquist频率为4KHz,而它的采样频率为5KHz,所以是满足Nyquist采样定律的,此时不会发生混叠。

运行效果如下:

Fs=1000sam/sec

这里使用的采样频率为1KHz,不满足Nyquist条件,因此会发生混叠。观察一下就会发生,1KHz采样得到的序列的频域波形和前面的频域波形不同,这就是混叠导致的,而且过低的采样率采集的信号的变换的不可逆的。

运行效果如下:



©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容

  • 一、傅立叶变换的由来 关于傅立叶变换,无论是书本还是在网上可以很容易找到关于傅立叶变换的描述,但是大都是些故弄玄虚...
    constant007阅读 4,493评论 1 10
  • 深入理解傅里叶变换Mar 12, 2017 这原本是我在知乎上对傅立叶变换、拉普拉斯变换、Z变换的联系?为什么要进...
    价值趋势技术派阅读 5,824评论 2 2
  • 姓名:张猛 【嵌牛导读】:如何对一个信号画出频谱并进行分析,从频谱中得到有用的信息 引用:http://blog....
    oliverabc阅读 18,047评论 0 6
  • 本文链接:个人站 | 简书 | CSDN版权声明:除特别声明外,本博客文章均采用 BY-NC-SA 许可协议。转载...
    虚胖一场阅读 9,514评论 0 6
  • 不要在疲乏烦躁慌乱痛苦的情况下做决定。
    鹿缤纷阅读 137评论 0 0