fft_my.m

来自「这是一个用于语音信号处理的工具箱」· M 代码 · 共 46 行

function x = fft_my(x,m,invrs)
      	
n=2^m;
y(1)=x(1);
for i=2:n
      	num=i-1;
      	j=0;
      	l=n;
      	for k=1:m
      		l=floor(l/2);
      		if num >= l 
      			j=j+2^(k-1);
      			num=num-l;
		end
      	end
      	y(i)=x(j+1);
end
pi=4.*atan(1.);
ldft=1;
ndft=n;
for k=1:m
        ldft=2*ldft;
        ndft=floor(ndft/2);
        for i=1:ndft
        	for j=1:floor(ldft/2)
        		arg=-(2.*pi*invrs)*(j-1.)/ldft;
        		w=cos(arg)+sqrt(-1)*sin(arg);
        		np=j+ldft*(i-1);
        		nq=np+ldft/2;
        		save=y(np)+w*y(nq);
        		y(nq)=y(np)-w*y(nq);
			y(np)=save;
		end
  	end	
end

for i=1:n
     	if invrs == 1 
		x(i)=y(i);
	end
	if invrs == -1 
		x(i)=y(i)/n;
	end
end