matlab 自己编写fft,基于matlab的自编fft算法以及可视化操作代码

该楼层疑似违规已被系统折叠 隐藏此楼查看此楼

f=input('请输入信号频率:')

N=input('请输入采样点数：')

T=input('请输入采样间隔:')

choicezero=input('增加零点? 请输入1/0 ')

n=0:T:(N-1)*T; %采样点

x=sin(2*pi*f*n);

if choicezero==1

a=input('请输入零点的个数!')

x=[x zeros(1,a)] %zeros 构建一个全为零的数组

N=N+a;

else

end %加零

k=0:N-1; %k取0-N-1共N个数

bianzhi=bi2de(fliplr(de2bi(k,length(de2bi(N-1)))))+1; %利用库函数进行变址运算，de2bi十进制转化为二进制，fliplr 矩阵的左右翻转，得到二进制数组

for l=1:N

X(l)=x(bianzhi(l)); %将采样后的值按照变址运算后的顺序放入X矩阵中

end

d1=1;

for m=1:log2(N)

d2=d1; %做蝶形运算的两个数之间的距离 d1=d1*2; %同一级之下蝶形结之间的距离

W=1; %蝶形运算系数的初始值

dw=exp(-j*pi/d2); %蝶形运算系数的增加量

for t=1:d2 %

for i=t:d1:N

i1=i+d2;

if(i1>N)break; %判断是否超出范围

else

p=X(i1)*W;

X(i1)=X(i)-p;

X(i)=X(i)+p; %蝶形运算

end

end

W=W*dw; %蝶形运算系数的变化

end

end

subplot(2,2,1);

t=0:0.0000001:N*T;

plot(t,sin(2*pi*f*t)); %画原曲线

subplot(2,2,2);

stem(k,x); %画采样后的离散信号

subplot(2,2,3);

stem(k,abs(X)/max(abs(X))); %画自己的fft的运算结果

subplot(2,2,4);

stem(k,abs(fft(x))/max(abs(fft(x)))); %调用系统的函数绘制结果