使用 MATLAB 进行离散傅立叶变换及其逆变换

随着 MATLAB 的出现及其带来的所有科学内置功能，复杂的工程场景发生了重大变化和简化。事实上，这种变化有助于帮助在科学领域（至少不是其他领域）从事技术研究的学生建立更好的可视化和实践技能。在这里，我们着眼于使用 MATLAB 实现一个基本的数学思想——离散傅立叶变换及其逆。

计算 DFT

定义离散傅立叶变换和逆变换如何将信号从时域转换到频域以及反之的标准方程如下：

DFT： $X(k)=\sum_{n=0}^{N-1}x(n).e^{-j2\pi kn/N}$对于 k=0, 1, 2....., N-1

国际电信联盟： $x(n)=\sum_{k=0}^{N-1}X(k).e^{j2\pi kn/N}$对于 n=0, 1, 2....., N-1

方程相当简单，可以简单地为求和执行重复/嵌套循环并完成它。但是，我们应该尝试使用另一种方法，即使用矩阵来找到问题的解决方案。许多读者会记得，时域/频域信号的 DFT 和 IDFT 可以用矢量格式表示如下：

$X(k) = \sum_{n=0}^{N-1}x(n).W_{N}^{kn}$

$x(n)=\frac{1}{N}\sum_{k=0}^{N-1}X(k).W_{N}^{-kn}$

$W_{N}=e^{-j2\pi/N}$

当我们将旋转因子作为矩阵的组成部分时，计算 DFT 和 IDFT 变得容易得多。因此，如果我们的频域信号是用X _N表示的单行矩阵，而时域信号也是用x _N表示的单行矩阵……

有了这个解释，我们需要做的就是创建两个数组，我们将在两个数组上发出矩阵乘法以获得输出。输出矩阵将始终是Nx1阶矩阵，因为我们将单行矩阵作为我们的输入信号（X _N或 x _N ）。这本质上是一个向量，为了方便起见，我们可以将其转置为水平矩阵。

算法（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

% 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

---

输出：

>> 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

% 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]

变换矩阵

时域序列