极客教程使用MATLAB进行离散傅里叶变换及其逆运算
随着MATLAB的出现以及它所带来的所有科学内置，复杂的工程场景发生了重大变化和简化。事实上，这种变化有助于为在科学领域从事技术研究的学生建立更好的可视化和实践技能，即使不是其他领域，也是如此。在这里，我们来看看如何使用MATLAB实现一个基本的数学思想–离散傅里叶变换及其逆向。
计算DFT
定义离散傅里叶变换和反转如何将信号从时域转换到频域，反之亦然的标准方程式如下。
DFT: 
for k=0, 1, 2….., N-1
IDFT: 
for n=0, 1, 2….., N-1
方程相当简单，人们可能只是简单地执行重复/嵌套循环的求和，然后就可以了。然而，我们应该尝试利用另一种方法，用矩阵来寻找问题的解决方案。许多读者会记得，时/频域信号的DFT和IDFT可以用矢量格式表示，如下所示。
当我们把孪生因子作为矩阵的组成部分时，计算DFT和IDFT就变得容易多了。因此，如果我们的频域信号是一个单行矩阵，用X N表示，时域信号也是一个单行矩阵，用X N表示 ……
有了这个解释，我们所需要做的就是创建两个数组，在这两个数组上进行矩阵乘法，以获得输出。由于我们将单行矩阵作为我们的输入信号（X N或x N），因此输出矩阵将始终是一个Nx1阶的矩阵。这本质上是一个矢量，为了方便，我们可以将其转置为一个水平矩阵。
算法(DFT)
1.获得输入序列和DFT序列的点数。
2.将获得的数据发送到一个计算DFT的函数中。声明一个新的函数并不是必须的，但代码的可读性和流程会变得更干净和明显。
3.使用length( )函数确定输入序列的长度，并检查长度是否大于点的数量。N必须总是等于或大于序列。如果你试图通过不满足条件来执行矩阵乘法，你会在命令窗口看到一个错误。
4.使用一个单独的数组对输入序列和N点的长度差异进行核算，该数组增加额外的零来拉长输入序列。这是用zeros(no_of_rows, no_of_columns)函数完成的，该函数创建了一个由零组成的二维数组。
5.根据输入的N值，创建W N矩阵。要做到这一点，需要实现2个 “for “循环–这是一个相当基本的程序。
6.简单地将已经创建的两个阵列相乘。这是一个所需频域信号样本的阵列。
7.通过内置函数abs(function_name)和angle(function_name)绘制出输出信号的幅度和相位。
% MATLAB code for DFT
clc;
xn=input('Input sequence: ');
N = input('Enter the number of points: ');
Xk=calcdft(xn,N);
disp('DFT X(k): ');
disp(Xk);
mgXk = abs(Xk);
phaseXk = angle(Xk);
k=0:N-1;
subplot (2,1,1);
stem(k,mgXk);
title ('DFT sequence: ');
xlabel('Frequency');
ylabel('Magnitude');
subplot(2,1,2);
stem(k,phaseXk);
title('Phase of the DFT sequence');
xlabel('Frequency');
ylabel('Phase');
function[Xk] = calcdft(xn,N)
L=length(xn);
if(N<L)
error('N must be greater than or equal to L!!')
end
x1=[xn, zeros(1,N-L)];
for k=0:1:N-1
for n=0:1:N-1
p=exp(-i*2*pi*n*k/N);
W(k+1,n+1)=p;
end
end
disp('Transformation matrix for DFT')
disp(W);
Xk=W*(x1.')
end
输出:
>> Input sequence: [1 4 9 16 25 36 49 64 81]
>> Enter the number of points: 9
算法(IDFT)
1.获得频域信号/序列作为输入（X（k））。这个序列的长度足以作为N（点）的值。
2.将这个数组传递给一个函数进行计算。
3.在函数中运行2个循环以创建矩阵。请注意，该矩阵在用于计算时必须是共轭的。你可以选择明确声明另一个数组，它是矩阵W N的共轭。
4.一旦创建了矩阵，用’*’获得共轭。 ‘并简单地与输入序列的转置相乘。我们需要转置，因为输入是一个行矩阵。当与我们创建的W N矩阵相乘时，W N的列数必须与X(k)的行数一致。
5.使用 stem(x_axis, y_axis) 绘制这个序列。不要使用plot( )，因为这不是一个CT信号。
% MATLAB code for IDFT
clc;
Xk = input('Input sequence X(k): ');
xn=calcidft(Xk);
N=length(xn);
disp('xn');
disp(xn);
n=0:N-1;
stem(n,xn);
xlabel('time');
ylabel('Amplitude');
function [xn] = calcidft(Xk) %function to calculate IDFT
N=length(Xk);
for k=0:1:N-1
for n=0:1:N-1
p=exp(i*2*pi*n*k/N);
IT(k+1,n+1)=p;
end
end
disp('Transformation Matrix for IDFT');
disp(IT);
xn = (IT*(Xk.'))/N;
end
输出:
>> Enter the input sequence: [1 2 3 4 5 9 8 7 6 5]
转化矩阵
时域序列