MATLAB Answers

Coding a function for the Fourier Transform.

3 views (last 30 days)
I have the following code for a function here,
function F = FourierTransInv(vec,arg)
% This function caluclates the Discrete Fourier Transform or the Inverse
% Fourier Transform.
% The input vec is a vector of length 2^n (n in N with 0) where the entries
% are all real numbers. arg is an input that determines whether or not the
% inverse or regular Fourier Transform is done. An arg of 1 corresponds to
% regular Fourier transform, while an arg of 0 is the inverse.
% The output F is the result of the Fourier or inverse transform.
if ~exist('vec')
error('Please make sure your inputs are valid.');
elseif ~isrow(vec) || ~isnumeric(vec)
error('Please make sure your inputs are valid.');
% Here 2^27 is the maximum because I observed that the difference in time
% between 2^20 and 2^21 is a factor of 2, where 2^20 took 30 seconds. So
% assuming the time increases by a factor of 2 each time, the minimum that
% 2^27 can take is 1 hour which is an extremely long time.
% This can of course be changed easily.
elseif ~ismember(length(vec),2.^(0:27))
error('Please make sure your input is a power of 2.');
elseif ~ismember(arg,[0,1])
arg=1;
warning('Argument is being set to 1 due to invalid input.');
end
N=length(vec);
F=ones(1,N);
switch arg
case 1
w=exp(-2*pi*1i/N);
case 0
w=exp(2*pi*1i/N);
end
switch N
case 1
F(1)=vec(1);
otherwise
Feven=FourierTransInv(vec(2:2:end),arg);
Fodd=FourierTransInv(vec(1:2:end),arg);
for k=1:N/2
F(k)=Fodd(k)+w^(k-1)*Feven(k);
F(k+N/2)=Fodd(k)-w^(k-1)*Feven(k);
end
end
if arg==0
% Division by sqrt(N) because this part always divides twice.
F=F/sqrt(N);
end
end
Now this works for example if I input the vector and then run it back through the inversion I get the same vector back. However, I tried some more complicated things such as the one produced with this line of code,
n=2^8; L=6*pi; dt=L/n; t=dt*[0:n-1];
y = sin(5*2*pi*t/L)+sin(10*2*pi*t/L);
I've tested many more than this and it all seems to lead to the same conclusion, if I run it through my function, taking the Fourier Transform it agrees with the function fft( . ) everywhere which is great. Then once I run it back through I get back almost the same thing except that every term is multiplied by a factor of , why is this? I cannot see anywhere in the code where this would happen, even less so because the simple works but not some slightly more complicated vector.

  0 Comments

Sign in to comment.

Accepted Answer

David Goodmanson
David Goodmanson on 22 Nov 2020
Edited: David Goodmanson on 22 Nov 2020
Hi Anthony,
the only problem here is the normalization. For the inverse transform you are dividing by sqrt(N) every time, (N being the 'local' N). However, the ifft has an overall factor of 1/2^N where 2^N is the length of the initial vector that is being transformed. So you need to divide by 2 at each level of recursion, but not if you have an element of length 2^0 = 1. If you replace
if arg==0
% Division by sqrt(N) because this part always divides twice.
F=F/sqrt(N);
end
with
if arg==0 & N~=1
F=F/2;
end
the results look good.

  0 Comments

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!