fft源代码(matlab)

对于正弦序列 )(8sin )(3n R n n x N ⋅⎪⎭⎫ ⎝⎛=π进行16点的dft 仿真

k=16;

n1=[0:1:15];

xa1=sin(2*pi*n1/k);

subplot(2,2,1)

plot(n1,xa1)

xlabel('t/T');ylabel('x(n)');

xk1=fft(xa1);xk1=abs(xk1);

subplot(2,2,2)

stem(n1,xk1)

xlabel('k');ylabel('X(k)');

n2=[0:1:7];

xa2=sin(2*pi*n2/k);

subplot(2,2,3)

plot(n2,xa2)

xlabel('t/T');ylabel('x(n)');

xk2=fft(xa2);xk2=abs(xk2);

subplot(2,2,4)

stem(n2,xk2)

xlabel('k');ylabel('X(k)');

function y=myfft(xr,n)

p=0:n-1 ; % 开始倒位序

nu=log2(n);

p1=p;

b=zeros(1,n);

for t=1:nu;

p2=floor(p1/2);

b=b*2+(p1-2*p2);

p1=p2;

end;

yr(p+1)=xr(b+1);

xr=yr; % 倒位序结束

t=0:n/2-1; %计算因子w 开始(只计算

for v=0:n/2-1;

w=exp(-2*i*pi*t/n);

end; %计算因子w 结束

for m=1:nu; % 计算x(k)开始

h=2^(m-1);

k=1;

while(k

for t=1:h;

y=bitshift(k-1,nu-m,nu)+1; %求w 的幂次数

xch(k)=xr(k)+w(y)*xr(k+h);

k=k+1;

end;

for t=1:h;

y=bitshift(k-1-h,nu-m,nu)+1; %求w 的幂次数

xch(k)=xr(k-h)-xr(k)*w(y);

k=k+1;

end;

end;

xr=xch;

end; % 计算x(k)结束

y=xr %输出变换后的结果