自建的离散傅里叶变换matlab程序实现及其与matlab自带函数比较举例

在matlab中有自带的离散傅里叶变换程序，即fft程序，但该程序是封装的，无法看到源码。为了比较清楚的了解matlab自带的实现过程，本文通过自建程序实现matlab程序，并与matlab自带的fft进行比较计算。

一、离散傅里叶变换的计算公式

在计算离散傅里叶变换的时候，通常会用到：
{ X ( k ) = ∑ n = 1 N [ x ( n ) ⋅ exp ⁡ ( − i ⋅ 2 π ( k − 1 ) ( n − 1 ) N ) ] s . t . { 1 ≤ k ≤ N } (1) \left\{ \begin{array}{l}X(k) = \sum\limits_{n = 1}^N {[x(n) \cdot \exp ( - i \cdot 2\pi \frac{{(k - 1)(n - 1)}}{N})]} \\s.t.\{ 1 \le k \le N\} \end{array} \right. \tag1 ⎩ ⎨ ⎧​X(k)=n=1∑N​[x(n)⋅exp(−i⋅2πN(k−1)(n−1)​)]s.t.{1≤k≤N}​(1)进行求解。
但有时会遇到所求解的向量长度N和变换过程中的长度K，大小不同。此时，会遇到$N\le K$的情况，和N>K的两种情况。

(1) 当$N\le K$时，则需要对向量$x$补零后，再离散傅里叶变换计算。计算公式为：

{ X ( k ) = ∑ n = 1 N p a d d e d [ x ( n ) ⋅ exp ⁡ ( − i ⋅ 2 π ( k − 1 ) ( n − 1 ) N p a d d e d ) ] s . t . { 1 ≤ k ≤ N p a d d e d } (2) \left\{ \begin{array}{l}X(k) = \sum\limits_{n = 1}^{{N_{padded}}} {[x(n) \cdot \exp ( - i \cdot 2\pi \frac{{(k - 1)(n - 1)}}{{{N_{padded}}}})]} \\s.t.\{ 1 \le k \le {N_{padded}}\} \end{array} \right. \tag2 ⎩ ⎨ ⎧​X(k)=n=1∑Npadded​​[x(n)⋅exp(−i⋅2πNpadded​(k−1)(n−1)​)]s.t.{1≤k≤Npadded​}​(2)
易知:$N\le N_{padded}=K$.

(2) 当N>K时，，则需要对向量$x$截断后，再离散傅里叶变换计算。计算公式为：

{ X ( k ) = ∑ n = 1 N t r u n c a t e d [ x ( n ) ⋅ exp ⁡ ( − i ⋅ 2 π ( k − 1 ) ( n − 1 ) N t r u n c a t e d ) ] s . t . { 1 ≤ k ≤ N t r u n c a t e d } (3) \left\{ \begin{array}{l}X(k) = \sum\limits_{n = 1}^{{N_{truncated}}} {[x(n) \cdot \exp ( - i \cdot 2\pi \frac{{(k - 1)(n - 1)}}{{{N_{truncated}}}})]} \\s.t.\{ 1 \le k \le {N_{truncated}}\} \end{array} \right. \tag3 ⎩ ⎨ ⎧​X(k)=n=1∑Ntruncated​​[x(n)⋅exp(−i⋅2πNtruncated​(k−1)(n−1)​)]s.t.{1≤k≤Ntruncated​}​(3)
易知:$N>N_{truncated}=K$.

二、基于上述理论编写myfft函数(matlab编程)

将自建的离散傅里叶变换的函数命名为myfft，编写程序如下：

function X=myfft(x,K)
% myfft函数根据傅里叶变换公式编写的离散傅里叶变换程序
% 输入
%     x：向量x
%     K: 变换后的向量X的长度

% 输出
%    X: 经过傅里叶变换得到的向量

% 变换依据：
% 对于长度为N的输入向量x，其离散傅里叶变换是长度为N的向量X,其具有元素：
%                    N
%      X(k) =       sum  x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
%                   n=1
%  myfft(x,K) 是一个K点的FFT，如果x小于K点，则补零后进行傅里叶变换；如果x大于K点，则截断后傅里叶变换。

%  by zddh and zsm
%  2023.10.24

N=length(x)

%% 1.如果x小于K点,补零运算
if N<=K
    x_padded=[x,zeros(1,K-N)];   %补零
    N_padded=length(x_padded);   %补零后的长度
    X=zeros(1,N_padded);         
for k=1:K
    for n=1:N_padded
    temp1=x_padded(n)*exp(-i*2*pi*(k-1)*(n-1)/N_padded);
    X(k)=X(k)+temp1;
    end
end
%% 2.如果x大于K点，则截断计算
else
   warning('K值小于N,则截断后进行傅里叶变换')
x_truncated=x(1:K);
N_truncated=length(x_truncated);
X=zeros(1,N_truncated)
for k=1:K
    for n=1:N_truncated
    temp2=x_truncated(n)*exp(-i*2*pi*(k-1)*(n-1)/N_truncated);
    X(k)=X(k)+temp2;
    end
end
end

三、自建的函数和matlab自带函数比较举例

(1) 编写程序

clc
clear all
close all
%% 1.构造将要变换的向量
dt=0.1
t=0:dt:10*pi;
x=sin(t)
N=length(x);
figure(1)
plot(t,x,'lineWidth',2)

%% 2.自建的离散傅里叶变换求解
K=200
X=myfft(x,K)
%% 3.matlab自带函数求解
X0=fft(x,K)


%% 4.比较
D_value=X-X0;

figure(2)
subplot(211)
plot(abs(X),'LineWidth',2)
hold on
plot(abs(X0),'LineWidth',2)
legend('myfft','matlabfft')
title('自建myfft和matlab自带函数fft比较')

subplot(212)
plot(abs(D_value),'LineWidth',2)
title('|X-X0|')

(2)运行结果

图1 生成的x向量

图2 使用两种方法结果
通过对图2两个子图观察比较可知，本文所编写的myfft函数和matlab自带的fft函数之间的误差非常小，在$10^{-12}$量级，同时验证了程序的理论公式(1)、（2）和（3）.