matlab离散傅里叶逆变换,手动实现离散傅里叶正变换与逆变换(程序+例子)

前言

实际野外采集到的信号都是有一定时间间隔的，因为本文主要关注于"离散傅里叶变换(DFT)"。另外，本文主要关注的是"如何使用"DFT这样一个工具(公式的使用+程序编写)，并非介绍其背后的原理。最后，本文将自己编写的程序与Matlab自带的相关语句的结果做对比，验证方法与程序的正确性。

DFT相关公式及使用

我们假设原始时域信号为x(n)，经DFT后频域信号为X(k)，其中n是采样点的序号，k是对应到频率域中的频率值，N是采样点总数。傅里叶正变换公式为：

公式1：DFT原始公式

注意1：在等号右边的最后一项

中，nk乘积是

的幂指数！

注意2：

是一个常数！其公式为：

公式1怎么用？先把它展开看看：

公式2：把公式1展开成多项来看

公式2看上去还是不够直观，我们再把它用矩阵表示成：

公式3：公式2的矩阵形式

是不是清晰多了呢？

已经知道

是一个常数！

那公式3"等号右边第一项"可进一步简化：

公式4：对公式3中右边第一项的简化

到此，DFT公式可变为：

这样的矩阵相乘。

有了上面的公式的基础，我们可以手写Matlab程序：

data = [2 3 8 9 12 4 8 10]; % 原始离散时域信号

N = length(data);

% DFT需要的相关矩阵:

WN = exp(-i*2*pi/N); % 常数

WN_nk = zeros(N)+WN; % WN_kn矩阵(初始)

xk = data'; % 时域信号振幅(列矩阵)

E = zeros(N); % 辅助的E(WN_kn的幂,单独拿出来算)

%%% 傅里叶正变换即结果 %%%

for row = 0:N-1

for cow = 0:N-1

E(row+1,cow+1) = row*cow;

WN_nk(row+1,cow+1) = WN_nk(row+1,cow+1)^E(row+1,cow+1);

end

end

Xk = WN_nk * xk; % 频域结果

fft(data); % Matlab自带语句,两者结果对比一样

IDFT相关公式及使用

很DFT的公式非常相似，只不过多了负号：

公式5：IDFT原始公式

将公式转换成矩阵，再写出矩阵E和对应的

即可。最后不要忘记除以总采样点数N。

对应的手写Matlab程序为：

% 继DFT代码后跟着写

% IDFT需要的相关矩阵：

WN_nk_n = zeros(N)+WN; % WN_kn矩阵(初始)

E_n = zeros(N); % 辅助的E(WN_kn的幂,单独拿出来算)

for row = 0:N-1

for cow = 0:N-1

E_n(row+1,cow+1) = -row*cow; % 注意负号

WN_nk_n(row+1,cow+1) = WN_nk_n(row+1,cow+1)^E_n(row+1,cow+1);

end

end

xk_n = real((WN_nk_n * Xk)/N)' % IDFT手写结果

ifft(Xk) % Matlab自带语句，上下俩结果对比一样。

更多完整程序与相关实例程序参考：傅里叶变换程序及实例.

补充说明：

频率域的优越性：时间域信号虽然可以反映原始信号不少的信息，但是这种图像不易修改！因为地球物理野外接收到的信号，不论是有效信号还是特定噪声，它们都各有自己所属的频带宽度！也就是说有效信号与噪声它们在频率域是很好识别和区分开来的！虽然频率信息也能少量反应在时域图像中(图像的紧密或舒张程度)，但是要想从时域图像中将二者分离开来是非常非常困难的！这也体现出频域的优越性。不只是信号提取与消除，频率域或者类频率域(后面的小波域、τp域等)还能做很多更高级的操作，这些都是时域图像所不能进行的。总结之：时域信号好采集不好处理，频域信号好处理但直接采集不了。