Reply by robert bristow-johnson●March 15, 20132013-03-15
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
Reply by clutchfft●March 15, 20132013-03-15
Does anyone have FFT implementation for Matlab? I'm not interested in the
built in Matlab function.