Prev: diffference scale adding lines in figure
Next: white gaussian noise
From: James Sampson on 3 May 2010 18:04
I wrote the following code to compute a FFT in MATLAB. I have it working, and have confirmed its accuracy using MATLAB's built in FFT function. My question is, how can I change this to do an inverse FFT?

CODE:
function f = myfft(p,data)
%1
if 2^(p+1)>length(data)
if floor(log(length(data))/log(2))==log(length(data))/log(2)
p=log(length(data))/log(2)-1;
else p=floor(log(length(data))/log(2));
end
fprintf('Due to an insufficient amount of data,\n')
fprintf('p has been changed to the maximum possible value: %d\n',p)
end
m=2^p;
M=m;
p=round(log(2*m)/log(2));
q=p-1;
zeta = exp(pi*i/m);
%2
C = data(1:2*m);
%3
xi=zeta.^(1:m);
xi=[xi -xi];
%4
K = 0;
%5
for L = 1 : p
%6
while K < 2*m-1
%7
for iter = 1 : M
%8
M1 = floor(K/exp(q*log(2)));
k = M1;
K2=0;
for count = 1 : p
K2 = 2*K2+k-2*floor(k/2);
k = floor(k/2);
end;
eta = C(K+M+1);
%9
if K2 ~= 0
eta = eta*xi(K2);
end;
C(K+M+1) = C(K+1)-eta;
C(K+1) = C(K+1)+eta;
%10
K = K+1;
end;
%11
K = K+M;
end;
%12
K = 0;
M = floor(M/2);
q = q-1;
end;
%13
while K < 2*m-1
%14
k = K;
sum=0;
for count = 1 : p
sum = 2*sum+k-2*floor(k/2);
k = floor(k/2);
end;
%15
if sum > K
temp = C(K+1);
C(K+1) = C(sum+1);
C(sum+1) = temp;
end;
%16
K = K+1;
end;
%17 and 18
C=conj(C);
ab=exp((-i*pi*m)*((1:2*m)-1)).*C(1:2*m)'/m;
a=real(ab(:));
b=imag(ab(:));
fprintf('Coefficients c(0), ... , c(2m-1)\n');
C
fprintf('\nCoefficients a(0), ..., a(m)\n');
a;
fprintf('\nCoefficients b(1), ..., b(m-1)\n');
b;
f=C;
%plotthisFFT(C,iter,begx,endx,steps)
end
function plotthisFFT(data,iter,begx,endx,steps)
y=zeros(1,steps);
x=begx:(endx-begx)/(steps-1):endx;
for iteration=1:steps
y(iteration)=(-4*A/pi)*sum(sin((1:iter).*pi/2).*(-1).^(1:iter).*(1:iter).^(-1).*(cos((1:iter).*pi/3)+1).*cos((1:iter).*2*pi*((x(iteration)-1/12))/T));
end
plot(x,y)
xlabel('t')
ylabel('y')
title('Fourier Series as given by ffseries')
end
 | 
Pages: 1
Prev: diffference scale adding lines in figure
Next: white gaussian noise