FFT函数，纯属个人学习兴趣，存在的问题，希望广大学者指正。
 
function Y = myfft(X,N)
%输入参数  X：输入信号
%          N： FFT的采样点数
%输出参数  Y：FFT结果，即频域信息

%判断信号的长度
if X<N   %若小于FFT的采样点数，则在末尾补0
    X = [X,zeros(1,N-length(X))];
else     %否则，截取前N个数作为FFT的输入
    X = X(1:N);
end

Y = zeros(1,N);  %初始化FFT结果的存储变量
for k = 1:N
    for n = 1:N
        y = X(n)*exp(-1i*2*pi*(k-1)*(n-1)/N);  %进行FFT转换
        Y(k) = Y(k)+y;    %将每次结果累加，得到最终的结果
    end
end
 
对应于上述FFT函数的IFFT函数
 
function X = myifft(Y,N)
%输入参数  Y：FFT结果，即频域信息
%          N： FFT的采样点数
%输出参数  X：IFFT后对应的时域信息

X = zeros(1,N);  %初始化IFFT结果的存储变量
for n = 1:N
    for k = 1:N
        x = real(1/N*Y(n)*exp(1i*2*pi*(k-1)*(n-1)/N));  %进行IFFT转换
        X(k) = X(k)+x;    %将每次结果累加，得到最终的结果
    end
end
 
验证：
 
clc
clear 

% 导入信号
load('EEG_signal_test.mat');  
signal = EEG_signal_test;  %原始信号
Fs = 160;    %采样频率
N = 256;  %FFT的采样点数

% 信号预处理
signal = signal(1:N)';
signal_nodc = signal-mean(signal);% 通过减去平均值来消除直流分量

%FFT变换
EEG_FFT = fft(signal_nodc,N);
EEG_FFT_Abs = abs(EEG_FFT);
EEG_myFFT = myfft(signal_nodc,N);
EEG_myFFT_Abs = abs(EEG_myFFT);
EEG_myFFT_Abs = EEG_myFFT_Abs';

%IFFT变换
EEG_IFFT = ifft(EEG_FFT,N);
EEG_myIFFT = myifft(EEG_myFFT,N);
EEG_myIFFT = EEG_myIFFT';

plot(signal);  hold on;  plot(EEG_IFFT)
 
验证结果截图
 
 ![在这里插入图片描述](https://img-blog.csdnimg.cn/2019101115365246.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3hxXzUyMA==,size_16,color_FFFFFF,t_7

 
上课和平时实验是用的FFT都是matlab自带的，最近需要在一个嵌入式处理器中对采集的信号进行FFT实时处理，所以想先在matlab中自己写一个fft，然后在底层去实现fft算法。
 
 
下面是matlab的代码（使用函数进行了封装）：
 
function [ ret_val ] = FFT1024( vector )
%UNTITLED8 此处显示有关此函数的摘要
%   此处显示详细说明
%UNTITLED7 此处显示有关此函数的摘要
%   此处显示详细说明
%======================================
%ret_val 为fft变换后返回的频域序列
%N 为点数
%vector 为变换前的序列
%======================================
vector_size = size(vector);
N = vector_size(2);
c = zeros(1,N);
%
%变址运算
%
j1 = 0;
for i = 1 : N
    if i < j1 + 1
        tmp = vector(j1 + 1);
        vector(j1 + 1) = vector(i);
        vector(i) =tmp;
    end
    k = N / 2;
    while k <= j1
        j1 = j1 - k;
        k = k / 2;
    end
    j1 = j1 + k;
end


%
%蝶形运算
%
%%%%%%%计算 N 的
dig = 0;
k = N;
while k > 1
    dig = dig + 1;
    k = k / 2;
end
%%%%%%
% m 为级; dist 为蝶形运两点的距离; n 为蝶形运算组数
%
n = N / 2;
for m = 1 : dig
    dist = 2 ^ (m - 1);
    idx = 1;
    for i = 1 : n
        idx1 = idx;
        for j1 = 1 : N / (2 * n)
            r = (idx - 1) * 2 ^ (dig - m);
            coef = exp(j * (-2 * pi * r / N));
            tmp                 = vector(idx);
            vector(idx)         = tmp + vector(idx + dist) * coef;
            vector(idx + dist)  = tmp - vector(idx + dist) * coef; 
            idx = idx + 1;
        end
        idx = idx1 + 2 * dist;
    end
    n = n / 2;
end


ret_val = vector;
end

 
测试代码：
 
clc;
clear;
AD_fre = 1023     %修改采样频率和采样点数，采样时间没有修改，通过修改改值可以观察周期截取对信号频谱的影响
Adc=2; %直流分量幅度
A1=3;   %频率F1信号的幅度
A2=1.5; %频率F2信号的幅度
F1=50; %信号1频率(Hz)
F2=75; %信号2频率(Hz)
Fs=AD_fre; %采样频率(Hz)
P1=0; %信号1相位(度)
P2=0; %信号相位(度)
N=AD_fre; %采样点数
t=[0:1/Fs:N/Fs]; %采样时刻
 
%信号
vector=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180);
%显示原始信号
figure(1)
plot(vector);
title('原始信号');
Y = FFT1024(vector); %做FFT变换
Ayy = (abs(Y)); %取模 
Ayy=Ayy/(N/2);   %换算成实际的幅度
Ayy(1)=Ayy(1)/2;
F=([1:N]-1)*Fs/N; %换算成实际的频率值，Fn=(n-1)*Fs/N
figure(2)
stem(F(1:N/2),Ayy(1:N/2));   %显示换算后的FFT模值结果
title('幅度-频率曲线图');
% 
% Pyy=[1:N/2];
% for i=1:N/2
% Pyy(i)=angle(Y(i)); %计算相位
% Pyy(i)=Pyy(i)*180/pi; %换算为角度
% end;
% subplot(414);stem(F(1:N/2),Pyy(1:N/2));   %显示相位图
% title('相位-频率曲线图');
 
 
 
 
 
 
 
 
 
 
 
 

matlab 快速傅里叶变换函数（fft)编写（横向读取数据进行变换）
 （欢迎指正）
 参考资料百度文库PPT：百度文库PPT
 
#matlab代码
 
#%快速傅里叶变换，不够2的整数幂的个数，末尾自动补齐0
function ret_val = myfft(Vector)
%因为输入的数据可能不是2的整数次幂，变换使得计算更加方便
[m,n]=size(Vector);%输入信号矩阵大小
num=ceil(log2(n));%向上取整
N=2^num;
vector=zeros(m,N);%申请足够大小矩阵
vector(:,1:n)=Vector(:,:);%将变换后的信号输入，不足补零
%变址运算
% temp=zeros(1,N);
% for p=1:m
%     for q=0:N-1
%         t=bin2dec(fliplr(dec2bin(q,num)));%变址运算
%         temp(1,q+1)=vector(p,t+1);
%     end
%     vector(p,:)=temp(1,:);
% end
%变址运算
for line=1:m%循环行数
    j1 = 0;
    for i = 1 : N%循环个数
        if i < j1 + 1
            tmp = vector(line,j1 + 1);
            vector(line,j1 + 1) = vector(line,i);
            vector(line,i) =tmp;
        end
        k = N / 2;
        while k <= j1
            j1 = j1 - k;
            k = k / 2;
        end
        j1 = j1 + k;
    end
end
for lines=1:m%循环行数
    for L=0:num-1%L表示运算等级或者层数
        dis=2^L;%dis表示奇偶组之间的距离
        for id=1:2^(num-L-1) %循环当前层数组数
            %进行同址运算
           for idx=1:dis%循环组内个数
               x1=(id-1)*2*dis+idx;%求得奇数数组的索引值
               x2=(id-1)*2*dis+dis+idx;%对应偶数数组的索引值
            temp1=vector(lines,x1)+vector(lines,x2)*W(L,(x1-1));%中间变量保存相应奇偶数组数据
            temp2=vector(lines,x1)-vector(lines,x2)*W(L,(x1-1));
            vector(lines,x1)=temp1;%所谓存入之前地址
            vector(lines,x2)=temp2;
           end
        end
    end
end
ret_val =vector;
function val=W(L,x)%旋转因子当层数为L,索引值为x
val=exp(-1j*2*pi*x/2^(L+1));
end
end
 
输入遥感影像数据进行测试结果截图*
 
 
输入矩阵或者数组进行测试并和matlab自带函数比较截图
 
 需要注意的是自己编写的myfft的横向读取数据和自动补零的操作方式和matlab自带fft函数有区别，不完全相同。

Matlab快速傅里叶变换程序（FFT）编写（不利用Matlab内置fft）
 
写一个基于基2FFT算法的傅里叶变换。利用原理：DIT-FFT（时域抽取法基2FFT）。
 DIT-FFT：
 设序列x(n)的长度为N，且满足N=2^M，M为自然数。
 ①按n的奇偶把x(n)分解为两个N/2的子序列。
 ②对x(n)做DFT，经过化简可以得到
 X(k)=X1(k)+W^k;
 X(k+N/2)=X1(k)+W^kX2(k);
 k=0,1,…,N/2-1;
 由此得到一个蝶形运算符号：
 
 ③这是DIT-FFT基本原理，按照N值的大小，还可进一步进行分解运算，可分解M次，将N点DFT分解为N个1点DFT和M级蝶形运算，而1点DFT就是时域序列本身。
 DFT与FFT运算量比较：
 N点DFT：复数乘法CM=NN
 复数加法CA=N(N-1)
 N点FFT： 复数乘法CM=MN/2
 复数加法CA=MN
 N=2^M
 代码：
 
function H=ditfft2(xn,M)
%DIT-FFT快速傅里叶变换程序
%程序名称：ditfft
%程序作者：grace_fight 2018/10/29

b=length(xn);                %调整补零
if(b<2^M)
    xn=[xn,zeros(1,(2^M)-b)];
end
b=length(xn);                %补零后xn长度
A=zeros(1,b);                %xn转换数组
N=2^M;                       %计算点数
nxd=bin2dec(fliplr(dec2bin([1:N]-1,M)))+1;%倒序排列序号
xn=xn(nxd);                               %倒序xn         
for i=1:N                                 %N个1点DFT，xn本身，赋值到数组A
   A(i)=xn(i);
end
for L = 1:M                               %DIT-FFT变换，M级蝶形变换
    B = 2^(L-1);                          %两个输入数据距离
    for J = 0:B-1;                        %旋转因子处理
        P=2^(M-L)*J;
        for k=(J+1):2^L:N;                %本次蝶形运算跨越时间
            W=exp(-j*2*pi*P/N);           %对应旋转因子
            T=A(k)+A(k+B)*W;              %进行蝶形运算
            A(k+B)=A(k)-A(k+B)*W;
            A(k)=T;
        end
    end
end
H=A;                                      %输出H为A
 
测试代码：
 
%测试程序
clear;                      %清屏
clc;
xn=[1 1 1 1];               %要转换的序列
M=5;                        %转换级数
B=fft(xn,32);               %matlab内置fft变换，用于对比
subplot(2,1,1);             %内置fft变换图像
[a1,b1]=size(B);
o1=1:b1;                    %内置fft变换坐标
stem(o1,abs(B),'.');        %绘图
H=ditfft2(xn,M);            %手工编写ditfft2函数
subplot(2,1,2);             %绘制手工编写的fft图像
[a1,b2]=size(H);
o2=1:b2;                    %ditfft2变换坐标
stem(o2,abs(H),'.')         %绘图
 
运行结果：
 
 参考：《数字信号处理第四版》高西全 丁玉美