function xk=fft_my(xn,N)
% *********************Declaration***************************
% File name: fft_my
% Author: @Harry
% Date: 2016/8/25 13:42:14
% Version Number: 1.0
% Modification history:[including time, version, author and abstract]
% 2016/8/25 version 1.0 xxx
% Abstract:
% 用时间抽取的基2FFT-DFT进行快速傅里叶变换
% *********************************end*******************************
%计算FFT的级数M
M=log2(N);
%对输入序列进行码位倒序
xn=bit_reverse(xn);
tempdata=xn;
for m=1:M
wcount=2^(m-1); %计算此级旋转因子的个数wcount
%计算这一级中wcount个旋转因子的值wvalue
for w=1:wcount
wvalue(w)=exp(-1i*pi*(w-1)/wcount);
end
%按照此级中旋转因子的个数分组
for num=1:wcount
group=(N/2)/(2^(m-1)); %同一旋转因子可分为group个组分别计算
for g=1:group
p=num+(g-1)*(2^m);
q=p+2^(m-1);
xk(p)=tempdata(p)+wvalue(num)*tempdata(q);
xk(q)=tempdata(p)-wvalue(num)*tempdata(q);
end
end
tempdata=xk;
end
function xk=bit_reverse(xn)
% *********************Declaration***************************
% File name: bit_reverse
% Author: @Harry
% Date: 2016/8/25 8:47:16
% Version Number: 1.0
% Modification history:[including time, version, author and abstract]
% 2016/8/25 version 1.0 xxx
% Abstract:
% 计算一个序列的码位倒序(序列长度为2^n)。把自然顺序的十进制转换成二进制数,
% 然后将这些二进制的首末尾倒序再重新转换成十进制,那么,这时的十进制
% 的排列就是码位倒序排列。(自然顺序必须从0开始)
% *********************************end***********************
%检查输入是否符合要求
[row,column]=size(xn);
if row~=1
error('输入必须是一维序列');
end
length=column; %序列长度
xk=zeros(1,length);
for i=0:length-1 %自然顺序必须从0开始
tempbinary=dec2bin(i,log2(length));
m=log2(length); %码位宽度
%开始对二进制的码位倒序
for j=1:m
tempdata(j)=tempbinary(m-j+1);
end
n=bin2dec(tempdata);
xk(i+1)=xn(n+1); %数组下标必须从1开始
end
快速傅里叶算法实现
最新推荐文章于 2024-03-01 16:23:09 发布