快速傅里叶变换_蝶形运算_按频率抽取基2-fft算法_MATLAB代码

function y=MyFFT_FB(x,n)

%MYFFT_TB:My Fast Fourier Transform Frequency Based

%按频率抽取基2-fft算法

%input:

% x -- 输入的一维样本

% n -- 变换长度，缺省时 n=length(x) 当n小于x数据长度时，x数据被截断到第n个数据

% 当n大于时，x数据在尾部补0直到 x 含n个数据

%output:

% y -- 1*n的向量，快速傅里叶变换结果

%variable define:

% N -- 一维数据x的长度

% xtem -- 临时储存x数据用

% m,M -- 对N进行分解 N=2^m*M,M为不能被2整除的整数

% two_m -- 2^m

% adr -- 变址，1*N的向量

% l -- 当前蝶形运算的级数

% W -- 长为 N/2的向量，记录 W(0,N),W(1,N),...W(N/2-1,N)

% d -- 蝶形运算两点间距离

% t -- 第l级蝶形运算含有的奇偶数组的个数

% mul -- 标量，乘数

% ind1,ind2 -- 标量，下标

% tem -- 标量，用于临时储存

%参考文献：

% 81c 输入参数个数检查

msg=nargchk(1,2,nargin);

error(msg);

%% 输入数据截断或加0

N=length(x);

if nargin==2

if N

xtem=x;

x=zeros(1,n);

x(1:N)=xtem;

N=n;

else % 截断

xtem=x;

x=xtem(1:n);

N=n;

end

end

%% 对N进行分解 N=2^m*M

[m,M]=factorize(N);

two_m=N/M;

%% 变换

if m~=0

%% 如果N可以被2整除

adr=address(m,M,two_m);

y=x; % 蝶形运算级数 l=m 时

%% 计算W向量

W=exp(-2*pi*i* ( 0:N/2-1 ) /N);

%% 蝶形运算

d=N/2;

t=1;

for l=1:m

% 加

for ii=0:t-1

ind1=ii*2*d+1;

ind2=ind1+d;

for r=0:d-1

tem=y(ind1)+y(ind2);

y(ind2)=y(ind1)-y(ind2);

y(ind1)=tem;

ind1=ind1+1;

ind2=ind2+1;

end

end

% 乘

for r=0:d-1

mul=W(r*t+1);

for ii=0:t-1

y(ii*2*d+d+1+r) = y(ii*2*d+d+1+r)*mul;

end

end

d=d/2;t=t*2;

end

%% 直接傅里叶变换

if M~=1 % N 分解含有非2因数M时，对y中每M个数据做直接傅里叶变换 for ii=1:two_m

y((ii-1)*M+1 : ii*M ) = DDFT( y((ii-1)*M+1 : ii*M ) );

end

%% 变址输出

y=y(adr+1);

else

%% 如果N 不能被2整除

y=DDFT(x);

end

end

%% 内嵌函数 ====================================================== function y=DDFT(x)

%% 直接离散傅里叶变换

%input:

% x -- 样本数据，N维向量

%output:

% y -- N维向量

%参考文献：

% 结构动力学，克拉夫，P82

% variable define

% s -- sum，用于求和

N=length(x);

y=zeros(size(x));

for n=1:N

s=0;

for m=1:N

s=s+x(m)*exp( -i*2*pi*(m-1)*(n-1)/N );

end

y(n)=s;

end

end

function [m,M]=factorize(N)

%% 对N分解

m=0;

while true

if mod(N,2)==0

m=m+1;

N=N/2;

else

break;

end

end

end

function adr=address(m,M,two_m)

%% 变址

% b -- 2^m * m 的矩阵，用来存储二进制数据

% ds -- 数，公差

adr=zeros(two_m,M);

b=de2bi(0:two_m-1,m);%转换为2进制注：matlab中二进制[0 1 1]=6 b=b(:,end:-1:1);% 逆序

adr(:,1)=bi2de(b);%2进制转换为10进制

if M~=1

ds=two_m;

adr=adr(:,1)*ones(1,M);

adr=adr+ds*ones(size(adr,1),1)*(0:M-1);

adr=reshape(adr',1,[]);

end

end