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

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

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

在计算离散傅里叶变换的时候，通常会用到：
$\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$ 进行求解。
但有时会遇到所求解的向量长度N和变换过程中的长度K，大小不同。此时，会遇到 $N\le K$ 的情况，和N>K的两种情况。

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

$\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$
易知: $\le {N_{padded}}=K$ .

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

$\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$
易知: $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|')