Caclulating the fourier transform
Mat Hunt
I am pretty used to matlab but have come across something which confuses me somewhat. I am performing a hydrodynamic stability analysis and I have written the wave profile as a fourier transform. I have calculated the functions I need to take the fourier transform of BUT I don't know how to go about doing it. As a test I tried to take the fourier transform of a gaussian which should give me a gaussian again.
The code I tried was:
>x=-3:0.001:3;
>y=exp(-x.^2);
f=abs(fft(y));
I get utter rubbish back, what am I doing wrong?
Cheers
Mat
Matt J
You're tripping over the circulant ordering properties of FFTs. You need to set up x so that the positive x is in the first half of your array and negative x is in the latter half. The following is a good way to do so:
N=6000; %number of samples
x= iffshift( -ceil((N-1)/2):N-1-ceil((N-1)/2) ) *.001;
y=exp(-x.^2);
f=fftshift(abs(fft(y)));
Mat Hunt
> "Mat Hunt" <hunt_...@hotmail.com> wrote in message <hf36si$f8...@fred.mathworks.com>...
Hi,
I typed in your code and received the following error message:
??? Undefined function or method 'iffshift' for input arguments of
type 'double'.
What am I doing wrong?
Matt J
> On Dec 1, 2:00?pm, "Matt J " <mattjacREM...@THISieee.spam> wrote:
> > "Mat Hunt" <hunt_...@hotmail.com> wrote in message <hf36si$f8...@fred.mathworks.com>...
> > > Hi,
> >
> >
> > > The code I tried was:
> > > >x=-3:0.001:3;
> > > >y=exp(-x.^2);
> > > f=abs(fft(y));
> >
> > > I get utter rubbish back, what am I doing wrong?
> >
> > You're tripping over the circulant ordering properties of FFTs. You need to set up x so that the positive x is in the first half of your array and negative x is in the latter half. The following is a good way to do so:
> >
> > N=6000; %number of samples
> > x= iffshift( -ceil((N-1)/2):N-1-ceil((N-1)/2) ) *.001;
> > y=exp(-x.^2);
> > f=fftshift(abs(fft(y)));
>
> Hi,
>
> I typed in your code and received the following error message:
> ??? Undefined function or method 'iffshift' for input arguments of
> type 'double'.
> What am I doing wrong?
A typo. It should be "ifftshift".
Since you've obviously not yet encountered the fftshift or ifftshift functions yet, you should type "help ifftshift" and "help fftshift" to see what they do.
Mat Hunt
I spotted the error. I went through the analysis and the answer was
still wrong. I get a peak at 3000, where I should get something less
than 5 centred at k=0.
When I take the fourier transform of a function f(x), I get the
function \hat{f}(\k). I want to be able to plot \hat{f}(\k) as a
function of k, is this so hard? I have a feeling that it shouldn't be
too hard and that matlab should be able to cope with it, otherwise it
will have to be an application of the trapizium rule....
Mat
Matt J
> Hi,
>
> I spotted the error. I went through the analysis and the answer was
> still wrong. I get a peak at 3000, where I should get something less
> than 5 centred at k=0.
That's because you've only told the plot() function what goes on the vertical axis, but not on the horizontal one. It therefore defaults to a horizontal axis of integers.
> When I take the fourier transform of a function f(x), I get the
> function \hat{f}(\k). I want to be able to plot \hat{f}(\k) as a
> function of k, is this so hard? I have a feeling that it shouldn't be
> too hard and that matlab should be able to cope with it, otherwise it
> will have to be an application of the trapizium rule....
I wouldn't say it's hard, but you have to know how the plot() function works and you have to know the subtle differences between an FFT and the continuous Fourier Transform. Below, I've done what I think gives you what you want. I've cranked up the number of samples N to get a smoother Gaussian curve:
N=1e6; %number of samples
dx=.001; %x axis sampling interval
df=1/N/dx; % f-axis sampling interval
%set up axes
Axis=( -ceil((N-1)/2):N-1-ceil((N-1)/2) );
xAxis=dx*Axis; %x-axis
fAxis=df*Axis; %frequency axis
%the actual data and transform operation
y=exp(-xAxis.^2); y=ifftshift(y);
f=abs(fft(y)); f=fftshift(f);
%%the plot
plot(fAxis, f)
xlim([-2,2])
xlabel 'Frequency (Hz)'
ylabel 'Fourier Transform (Approx)'
Matt J
> Hi,
>
Another thing you should realize is that different professions do not agree on what the formula for the Fourier transform should be, and in particular on what scaling constants it should have. So whether you would expect something less than 5 at k=0 depends on whether you are a physicist, engineer, or whatever.
According to my profession, matching the scale of the FFT to the continuous Fourier transform requires the following scaling:
f=abs(fft(y))*dx;
That is, to approximate the continuous Fourier transform (an integral) by the FFT (a sum), you must multiply the sum with the x-axis sampling interval dx. However, if you are a physicist, additional scale factors might be necessary to get what you want.
The best thing is to look at "help fft" which shows the formula for FFT used by MATLAB. Then, you can analyze how it should be scaled to meet your tastes.
Mat Hunt
I have gotten the scaling as I want them, I will have to keep tract of all these 2\pi's floating around everywhere. It's going to be a pain but not too much of a hardship.
Mat
Mat Hunt
Mat
Matt J
> Turns out that I require the inverse Fourier transform, so I guess it is pretty much the same but replace fft with ifft?
>
And you'll need to divide by dx instead of multiplying by it.
Mat Hunt
I modified your script to calculate the inverse fourier transform and it worked wonderfully, I then went on to calculate the inverse Fourier transform that gives me the shape of my free surface and I ran into problems, wondered if you could tell me where I am making my errors.
I tried to use functions which were bounded (sech and tanh) to makes things easier for matlab but all I get back is gibberish and I don't know why.
The script I used to compute the inverse Fourier transform is:
%This sees if I can getthe inverse Fourier transform working in good old
%matlab.
N=1e6; %number of samples
dx=.001; %x axis sampling interval
df=1/N/dx; % f-axis sampling interval
%set up axes
Axis=( -ceil((N-1)/2):N-1-ceil((N-1)/2) );
xAxis=dx*Axis; %x-axis
fAxis=df*Axis; %frequency axis
%the actual data and transform operation
y=sech(pi^2*xAxis); y=ifftshift(y);
%Now try and reverse the process and calculate the inverse of the fourier
%transform.
F=pi*sech(pi^2*fAxis);
F=ifftshift(F);
q=abs(ifft(F));
q=fftshift(q)/dx;
plot(xAxis,q,'r');
axis([-2 2 0 1]);
hold on
pause(2)
plot(xAxis,sech(xAxis));
The script that I used to find the shape of my free surface is:
%This plots the free surface of the wave.
B=1;
E_b=1;
U=1;
h=1;
N=1e6; %number of samples
dx=.001; %x axis sampling interval
df=1/N/dx; % f-axis sampling interval
%set up axes
Axis=( -ceil((N-1)/2):N-1-ceil((N-1)/2) );
xAxis=dx*Axis; %x-axis
fAxis=df*Axis; %frequency axis
%the actual data and transform operation
x=sech(pi^2*fAxis).^2;
y=tanh(fAxis*h);
w=(B/U)-(E_b*fAxis/U)+(4*pi^2*fAxis.^2/U);
v=-4*pi^2*U*fAxis+w.*y;
A=-2*pi^2*i*x./v;
%Now try and reverse the process and calculate the inverse of the fourier
%transform.
F=-i*A.*sinh(fAxis*h)/(2*pi*U);
F=ifftshift(F);
q=abs(ifft(F));
q=fftshift(q)/dx;
plot(xAxis,q,'r');
Can you see what I am doing wrong? Are my numbers too big?
Matt J
> Can you see what I am doing wrong? Are my numbers too big?
Both sets of code give me nice and very similar looking curves. I can't see anything wrong, but then again, I don't know what I should be seeing :)
Mat Hunt
Gut feeling is that the two curves shouldn't be the same. The answer should show the response of a free surface to a pressure distribution p(x)=sech(x). So I don't think that the response should be the same. I think (from a research paper that my supervisor wrote) that I should see a series of beats.
I will have a go at doing the transform numerically and see if I get the same answer.
Mat
Matt J
> Odd, I got the first program to work quite easily but I kept getting a vector of NaN's for the second program. This was the reason why I posted the second program
>
> Gut feeling is that the two curves shouldn't be the same. The answer should show the response of a free surface to a pressure distribution p(x)=sech(x). So I don't think that the response should be the same. I think (from a research paper that my supervisor wrote) that I should see a series of beats.
>
> I will have a go at doing the transform numerically and see if I get the same answer.
I think I found the problem. In this computation
A=-2*pi^2*i*x./v;
one component of v is zero. The division by 0 gives a non-finite value and propagates into the ifft as NaNs.
Jorian
Thanks Matt!
Best regards, Jorian Seokaner!
http://www.mathworks.com/matlabcentral/newsreader/author/126323