基-2fft算法的软件实现 matlab代码

% 基于Matlab的时间抽取基2FFT算法

function y=myditfft(x)

%本程序对输入序列实现DIT-FFT基2算法，点数取大于等于长度的2的幂次

%------------------------------------

% Leo's fft program(改编网上的一个程序)

%------------------------------------

m=log2(2^nextpow2(length(x))); %求的x长度对应的2的最低幂次m

N=2^m;

if length(x)<N

x=[x,zeros(1,N-length(x))]; %若长度不是2的幂，补0到2的整数幂

end

x;

%--------------------------------------------------------------------------

%对输入序列进行倒序

%如果输入序列的自然顺序号I用二进制数（例如n2n1n0）表示

%则其倒位序J对应的二进制数就是（n0n1n2），这样，在原来自然顺序时应该放x(I)的

%单元，现在倒位序后应放x(J)。

%--------------------------------------------------------------------------

%以下程序相当于以下程序:

%nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1; %求1：2^m数列的倒序

%y=x(nxd); %将倒序排列作为初始值

%--------------------------------------------------------------------------

NV2=N/2;

NM1=N-1;

I=0;

J=0;

while I<NM1

if I<J

T=x(J+1);

x(J+1)=x(I+1);

x(I+1)=T;

end

K=NV2;

while K<=J

J=J-K;

K=K/2;

end

J=J+K;

I=I+1;

end

x;

%--------------------------------------------------------------------------

%以下程序解释:

%第一级从x(0)开始,跨接一阶蝶形,再取每条对称

%第二级从x(0)开始,跨接两阶蝶形,再取每条对称

%第m级从x(0)开始,跨接2^(m-1)阶蝶形,再取每条对称

%--------------------------------------------------------------------------

for mm=1:m %将DFT做m次基2分解，从左到右，对每次分解作DFT运算

Nmr=2^mm;

u=1; %旋转因子u初始化

WN=exp(-j2pi/Nmr); %本次分解的基本DFT因子WN＝exp(-i2pi/Nmr)

for n=1:Nmr/2 %本次跨越间隔内的各次碟形运算

for k=n:Nmr:N %本次碟形运算的跨越间隔为Nmr=2^mm

kp=k+Nmr/2; %确定碟形运算的对应单元下标(对称性)

t=x(kp)u; %碟形运算的乘积项

x(kp)=x(k)-t; %碟形运算的加法项

x(k)=x(k)+t;

end

u=uWN; %修改旋转因子，多乘一个基本DFT因子WN

end

end

y=x; %输出

N=6000;

fs=12000000;

n=0:N-1;

%%%num为6000000个数据值;

F=fft(num);

f=n/Nfs;

plot(f,abs(F));

fs=200;

N=1024;

n=-N:N-1;

t=1/fsn;

f=fs/Nn/2;

x=exp(-t)(t>=0)-exp(t)(t<0);

y=fftshift(fft(x,2N));

plot(f,abs(y));

fft结果出来的是从0频到fs的频谱

要画成从-fs/2到fs/2的频谱，需要再用fftshift函数

貌似这是个书呆子编滴程序，以为乘法有交换律就随意改变位置，完全不管物理意义好不好理解，其实不就是采样定理要除2嘛！fft后最高频率是采样频率滴一半，计算频率数组时是以此原理再按几分之几均分其它部分，其实最易理解滴是(0:length(Yc)-1) /length(Yc)(Fs/2)这种表达，只不过原式与这式数学上计算结果相同。还有乘除法不用点运算，小心编译出错，真是相当糟糕的代码！

如果没看错是再求信号的功率谱，

Fs=2048;

%

设定采样点数

FFTN=1024;

%设定FFT点数

t=(0:1023)/Fs;

%

设定正弦信号的采样时间点

x=2sin(2pi256t+pi/8);

%

产生正弦信号

y=fftshift(fft(x,FFTN));

%这句没太看懂，应该是求信号的功率谱，但结果是复数，需取模，这句需楼主在核实下

yam=abs(y)2/FFTN;

%需取模，

f=((0:FFTN-1)-FFTN/2)Fs/FFTN;

%设定画图时的横坐标，Hz为单位

figure(1);plot(f,yam);xlabel('f(Hz)');ylabel('幅值');

%

以横坐标Hz，纵坐标取模结果，画图，

[amax

index]=max(yam);

%取纵坐标最大值和序号

gmax=angle(y(index));

%这句也没太看，取出模值最大值进行什么处理，

fmax=f(index);

%取出模值最大点对应的频率值

以上除了标注的，其他应该没问题，楼主在看看那两句吧，查查matlab的help

现成的FFT程序网上可以下到很多，但如何定义相似度需要考虑一下。可以考虑用相对变化率的形式（相对变化率=（数值1-数值2）/数值1 或者 相对变化率=（数值1-数值2）/数值2 ）。

有一种情况跟你的需求很像：设计滤波器后说明滤波效果。这种情况，需要比较变化前后的信号进行幅频特性、相频特性曲线，以此说明滤波效果。在这一过程中就需要对两个信号分别进行FFT变换，以求得幅频曲线和相频曲线。

具体过程如下：

step1：将横坐标定义为t（matlab赋值语句t=[,,];），将纵坐标定义为y（matlab赋值语句y=[,,];），采样频率就是临近两个横坐标差值的倒数（一般横坐标都应为时间）；

step2：新建一个m文件（快捷键为ctrl+N）；

step3：将如下程序复制到m文件中

t=[,,];%填入横轴数据

y=[,,];%填入纵轴数据

N = size(t,2);%行向量时，列向量是为N = size(t,1);

N = 2^(nextpow2(N)-1);

Y = fft(y,N);

mag_Y = abs(Y)/N2;%各个频率点处的幅值，这个就是fft变换后的数据

f = fs/2linspace(0,1,N/2+1);%对应的频率值

plot(f,mag_Y(1:N/2+1))%显示变换后的曲线

step4：亲测实例

t = 1:0001:10;

y= 2sin(2pi100t);%频率为100Hz

%=================以下同step3中代码==============

  N = size(t,2);%行向量时，列向量是为N = size(t,1);

N = 2^(nextpow2(N)-1);

Y = fft(y,N);

mag_Y = abs(Y)/N2;%各个频率点处的幅值，这个就是fft变换后的数据

f = fs/2linspace(0,1,N/2+1);%对应的频率值

plot(f,mag_Y(1:N/2+1))%显示变换后的曲线 

%=================以上同step3中代码==============

结果：

由于matlab自带FFT函数有些小瑕疵，所以在100Hz处幅值虽然接近2，但还是有些偏差的，对于偏差的修正网上也有相应的方法，如果需要在留言。

用MATLAB 实现傅里叶变换:

用户任意输入一个函数，然后，输出函数的傅里叶变换函数，然后输出振幅频率 。

x=sin(2pit); %任意输入一个函数。

y=fft(x); %傅里叶变换函数。

plot(abs(y)); %振幅频率。

以上就是关于基-2fft算法的软件实现 matlab代码全部的内容，包括:基-2fft算法的软件实现 matlab代码、在MATLAB中怎么对数据做fft、matlab fft 用法等相关内容解答，如果想了解更多相关内容，可以关注我们，你们的支持是我们更新的动力！