Matlab FFT Code

Does anyone have FFT implementation for Matlab?  I'm not interested in the
built in Matlab function.  

On 3/15/13 11:15 AM, clutchfft wrote:
> Does anyone have FFT implementation for Matlab?  I'm not interested in the
> built in Matlab function.

since your e-dress looks like crap, i guess i'll just append it below. 
it short enough.  (i wrote this for a demonstration about 5 years ago of 
what happens to the FFT when the quantization word gets short.  the 
files that put in the quantization are other files.)  this code was not 
meant to be fast, but just to be able to show what happens in each 
butterfly, group, and pass.

you will have to worry about unwrapping long lines, since i cannot turn 
off the word wrap in my email client.  as far as i know, it worked the 
last time i used it.  lemme know if it doesn't.


-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."


FILE: COEF_TABLE.M

function coefs = coef_table(x)

	half_N = length(x)/2;
	
	if round(log2(half_N)) == log2(half_N)
	
		coefs = exp(i*linspace(0, (half_N-1)/half_N*pi, half_N));

	else

		error('length of argument must be power of two');
	
	end



FILE: FORWARD_FFT.M

function y = forward_fft(x, coefs)

%
%	length(x) must be an integer power of 2.
%	coefs = coef_table(x) must be called first.
%	input x is in normal order, output y is in bit-reversed order.
%	call y = bit_reverse(y) afterward to swap output back to normal order.
%

	y = x;
	
	N = length(x);
	p = log2(N);
	
	if round(p) == p
		
		bfly_length = N/2;
		n_group = 1;
	
		for pp = 1:p
			
			group_top = 1;
			for g = 1:n_group
				
				top = group_top;
				bottom = top + bfly_length;
				coef_ptr = 1;
				for b = 1:bfly_length								% num bflys = bfly_length
					
					A = y(top) + y(bottom);
					B = coefs(coef_ptr)*(y(top) - y(bottom));
					
					y(top) = A;
					y(bottom) = B;
					
					top = top + 1;
					bottom = bottom + 1;
					coef_ptr = coef_ptr + n_group;
				end
			
				group_top = group_top + 2*bfly_length;
			end
			
			bfly_length = 0.5*bfly_length;
			n_group = 2*n_group;
			
		end
		
	else

		error('length of argument must be power of two');
	
	end



FILE: BIT_REVERSE.M

function y = bit_reverse(x)

	y = x;
	
	N = length(x);
	p = log2(N);
	
	if round(p) == p
	
		for n = 0:(N-1)
			r = 0;
			m = n;
			for pp = 1:p
				r = bitshift(r, 1);
				r = bitor(r, bitand(m,1));		% reverse the bits
				m = bitshift(m, -1);
			end
		
			y(r+1) = x(n+1);	
	
		end
	
	else

		error('length of argument must be power of two');
	
	end



FILE: INVERSE_FFT.M

function y = inverse_fft(x, coefs)

%
%	length(x) must be an integer power of 2.
%	coefs = coef_table(x) must be called first.
%	input x is in bit-reversed order, output y is in normal order.
%	call x = bit_reverse(x) first to swap input to bit-reversed order.
%

	y = x;
	
	N = length(x);
	p = log2(N);
	
	if round(p) == p
		
		coefs = conj(coefs);
		
		bfly_length = 1;
		n_group = N/2;
	
		for pp = 1:p
			
			group_top = 1;
			for g = 1:n_group
				
				top = group_top;
				bottom = top + bfly_length;
				coef_ptr = 1;
				for b = 1:bfly_length								% num bflys = bfly_length
					
					A = 0.5*y(top);
					B = 0.5*y(bottom)*coefs(coef_ptr);
					
					A = A + B;
					B = A - 2*B;

					y(top) = A;
					y(bottom) = B;
										
					top = top + 1;
					bottom = bottom + 1;
					coef_ptr = coef_ptr + n_group;
				end
			
				group_top = group_top + 2*bfly_length;
			end
			
			bfly_length = 2*bfly_length;
			n_group = 0.5*n_group;
			
		end
	
		coefs = conj(coefs);
		
	else

		error('length of argument must be power of two');
	
	end