FFTW questions

[Biftheunderstudy]

I'm working with Intel's FFTW libraries and I have a few questions about phases.

First off, as a test example I'm taking the fft of a gaussian pulse. I've set it up so that the pulse is centered at x=5, with x ranging from 0 to 10.

I plot the initial pulse, so far so good, its the gaussian I expected.
Next, I take the fftw. Now I get the phase shifted version (since the dft presumes that the first array value is the origin).

Finally, I fft back and happily get my initial gaussian back (after dividing by n to normalize).

Now there are a couple of complications...
I wrote a function to shift the data back to the center after a fft, this means that when I do the forward transform and plot the results, I get a gaussian again instead of the phase shifted one. (In practice this means I take the second half of the array and set it to the beginning of the array, there are complications if the array is odd or even). Anyone have suggestions for how to write a good shifting routine that is symmetric for odd data?

So far so good, now I get a nice gaussian for all three plots and when I actually go to put this in my code I can work on the transformed data.

Now here is the problem, if the pulse is not centered, say its at x=3. When I forward transform I get the usual, then I backwards transform the pulse is now centered at x = 8! (Same thing happens if I shift in between the forward and backward transforms)

Anyone have any advice on how to take the fft so that I can take the forward transform, do some work, then take the backwards transform and not have to deal with phases messing things up?

 

[TecHNooB]

what do you mean by shift the data back to the center after a fft?

when you take the fft, you interpret what you get as [0 to fs/2, -fs/2 to 0] repeating.

when you take the FFT of a gaussian, you get a gaussian. in fact, you get a repeating gaussian. but it's not represented as -fs/2 to 0 to fs/2. it's shifted by fs/2. So you're going to see the right half of the gaussian first.. then the left half. It should look like a bathtub. a phase shifted gaussian would still give you the bathtub, but your phase components will be shifted by however much the gaussian is shifted in time. the magnitude spectrum will still look like a gaussian about zero (or a bathtub if you didn't shift by fs/2, which is the default).

I thought this was familiar.
http://forums.anandtech.com/showthread.php?t=2128793&highlight=

 

[Biftheunderstudy]

Yes, I understand that. I'm now using the Intel fast fourier transform libraries in my own code. I didn't think it would be proper to dig up an old thread when this one has a different question.

So, if I'm understanding this correctly then:
If I take my gaussian centered at x=3 (x=0 to x=10) and take the forward transform, I'll get the bathtub shape. I can then shift this once to get the gaussian shape back, centered at x=3 (or rather k=3). Then shift back to the gaussian shape again, and do the backwards transform. Now however, the gaussian is centered on x=8, even if I skip the shifting.

Does this mean that the data must be centered before I can do any fft's? This might be useful information because, in general, I won't be dealing with symmetric data.

 

[TecHNooB]

Yes, I understand that. I'm now using the Intel fast fourier transform libraries in my own code. I didn't think it would be proper to dig up an old thread when this one has a different question.

So, if I'm understanding this correctly then:
If I take my gaussian centered at x=3 (x=0 to x=10) and take the forward transform, I'll get the bathtub shape. I can then shift this once to get the gaussian shape back, centered at x=3 (or rather k=3). Then shift back to the gaussian shape again, and do the backwards transform. Now however, the gaussian is centered on x=8, even if I skip the shifting.

Does this mean that the data must be centered before I can do any fft's? This might be useful information because, in general, I won't be dealing with symmetric data.

you can leave the data alone. just learn how to interpret the indicies of the fft.

normally ppl shift by fs/2 so they can interpret the indicies as -fs/2 to fs/2. your gaussian would be in the middle bin, not x = 3. there are some nuances depending on whether the FFT is for an even or odd number of points.

i'll write something in matlab for you, don't run away. ignore the frequency axis for bathtub, they're wrong -.- hm.. some of my comments are hella bad too. sample frequency should be sample period, but w/e.

clc, clear all, close all

% Sampling stuff.
N = 300; % # of samples.
x_start = -10;
x_end = 20;
X = x_end - x_start; 
dx = X / (N - 1); % Sample frequency.
x = x_start:dx:x_end;

% Frequency
fs = 1/dx;

if(mod(N,2) == 0) % Even number of samples.
    df = fs / N;
    f = -fs/2:df:fs/2-df;
else % Odd number of samples.
    df = fs / (N - 1);
    f = -fs/2:df:fs/2;
end

% Normal gaussian.
mean = 0;
sigma = 1;
g = 1/sqrt(2*pi*sigma^2)*exp(-(x - mean).^2/2/sigma^2);

figure(1);
subplot(4,1,1);
stem(x,g);

G = fft(g);
subplot(4,1,2);
stem(abs(G));
title('Bathtub');
xlabel('Frequency');
ylabel('Magnitude');
axis tight;

subplot(4,1,3);
stem(f,fftshift(abs(G)));
title('Gaussian');
xlabel('Frequency');
ylabel('Magnitude');
axis tight;

subplot(4,1,4);
stem(f,fftshift(angle(G)));
title('Gaussian');
xlabel('Frequency');
ylabel('Phase');
axis tight;

% Shifted gaussian
for mean = linspace(-5,15,100);
    sigma = 1;
    g = 1/sqrt(2*pi*sigma^2)*exp(-(x - mean).^2/2/sigma^2);

    figure(2);
    subplot(4,1,1);
    plot(x,g,'.');

    G = fft(g);
    subplot(4,1,2);
    plot(abs(G),'.');
    title('Bathtub');
    xlabel('Frequency');
    ylabel('Magnitude');
    axis tight;

    subplot(4,1,3);
    plot(f,fftshift(abs(G)),'.');
    title('Gaussian');
    xlabel('Frequency');
    ylabel('Magnitude');
    axis tight;

    subplot(4,1,4);
    plot(f,fftshift(angle(G)),'.');
    title('Gaussian');
    xlabel('Frequency');
    ylabel('Phase');
    axis tight;
    
    pause(0.01);
end

 

[TecHNooB]

you know what, i just re-read your op and realized i didnt answer your question and repeated half of what you said. what filter did you apply to the gaussian? you probably phase shifted it unintentionally. if whatever frequency operation you're doing shifts your pulse then that's what it does, nothing is necessarily wrong. if you dont want the phase, just time shift it back yourself or add that into the filter's frequency response.

man.. i need to sleep.

 

[Biftheunderstudy]

Maybe a little background will help clarify.

I've been solving some DE's and I thought it might be a good idea to Fourier transform the lot of them so that I can some them more effectively. I've done most of the work analytically and I'm solving the resulting equations numerically.

Somewhere along the line, I'll need to transform the numeric solution back to real space to look at the spatial distribution of one of the intermediate steps.

Also, I've fft'ed my initial conditions which happen to be a gaussian. I'm using the Intel FFTW library, so it does the most straight forward discrete transform with a bunch of optimizations, so no filters or anything applied. This normally produces the "bathtub" so I shifted it before letting the rest of the code at it.

My worry is that when the solver does its thing and I want to look at the output spatially that it wont be centered or symmetric and the shift + backwards fft will put details in the wrong orientation. So, I would like to be able to nail down how this works in the general asymmetric case without the requirement that the whole thing be centered.

Thanks for the help by the way, no one that I work with deals with discrete fft's so I'm on my own here.

 

[TecHNooB]

nm. can you post more detail?

are you doing something like this? http://www.math.ubc.ca/~feldman/m267/pdeft.pdf

i've never actually used fourier transforms to solve differential equations numerically. however, one way that might work is if you have some transfer function H(x,f) and initial conditions a(x,0); you can try to take the inverse of H(x,f) straight away and simply convolve with your initial conditions for all x. convolution is usually a good way to avoid symmetry issues cuz its done in the time domain.

 

[esun]

Why is symmetry even relevant here? If you have some signal x[n] and you compute X(jw), then inverse transform to get x[n] again, it does not matter if x[n] is symmetric to begin with. If you're getting some shift at the end, then you're doing something wrong when you're shifting indices in the frequency domain.

 

[TecHNooB]

o man i really needed sleep last night. even my post about needing sleep was full of wrongness. i see what you're saying. so what you did was you took the fft of a gaussian centered at x = 3 (for x = 0 to 10) and you got the bathtub. that part is fine. but when you shift it, you shifted it to k = 3? the gaussian should always be centered in the middle (k = 5 in your case). if it's not, you gotta interpret your frequency indexes very carefully. so yes, do see my example matlab code

you need to pay close attention to the phase information, not the magnitude information. that's what's conveying the time shift (or space shift, w/e). the magnitude info is going to stay roughly the same.

the first index is 0 hz or DC or average value or w/e. knowing that, you can shift however you want but when you multiply other frequency responses together, you have to match frequencies. the inverse fft assumes you're going to give it the 0 hz to fs/2, -fs/2 to 0 hz frequency response, so after whatever processing you do, you have to shift the response around such that you give the inverse fft the response with the frequency indexes it expects.

 

[Biftheunderstudy]

TechNoob, I did take a look at the matlab script, it did help in the sense that it at least proves that I'm implementing the libraries correctly and gave a little insight into how the individual transforms are storing the data.

I've sort of hacked together a subroutine that makes things work. What I did was took the fft of my data/signal/array, then when I want to inverse transform it I:

ifft, then fft again, then ifft again and you get back what you started with in the proper spatial positioning.

esun, what I'm saying is:
If I have a gaussian, it's centered on x=3, where x ranges from 0 to 10.
Now I take an fft of it using either matlab's fft function or the fftw c libraries or what ever, then I take the inverse transform, ifft in matlab or just switching the fftw_forward to fftw_backward in the libraries. You would expect that you should get your same gaussian back that you started with; you don't. You get a gaussian that is now centered at x=7.

My hunch is that if you took a gaussian centered in the middle and labeled the left side red and the right side blue that after a fft->ifft you would end up with the reverse situation (red on the right, blue on the left).

To fix this, I just did an extra set of fft->ifft to get everything back the way I expect them. To reiterate, I'm not doing signal or image processing, I'm solving ode's using spectral methods so symmetry might matter. I'm not actually sure what the effects of having time shifted signals on the solutions to these equations might be, I would rather not have to worry about it in the first place.

In matlab code this would be something like ifft(fft(ifft(fft(g)))) to get the original 'g' back.

 

[TecHNooB]

TechNoob, I did take a look at the matlab script, it did help in the sense that it at least proves that I'm implementing the libraries correctly and gave a little insight into how the individual transforms are storing the data.

I've sort of hacked together a subroutine that makes things work. What I did was took the fft of my data/signal/array, then when I want to inverse transform it I:

ifft, then fft again, then ifft again and you get back what you started with in the proper spatial positioning.

esun, what I'm saying is:
If I have a gaussian, it's centered on x=3, where x ranges from 0 to 10.
Now I take an fft of it using either matlab's fft function or the fftw c libraries or what ever, then I take the inverse transform, ifft in matlab or just switching the fftw_forward to fftw_backward in the libraries. You would expect that you should get your same gaussian back that you started with; you don't. You get a gaussian that is now centered at x=7.

My hunch is that if you took a gaussian centered in the middle and labeled the left side red and the right side blue that after a fft->ifft you would end up with the reverse situation (red on the right, blue on the left).

To fix this, I just did an extra set of fft->ifft to get everything back the way I expect them. To reiterate, I'm not doing signal or image processing, I'm solving ode's using spectral methods so symmetry might matter. I'm not actually sure what the effects of having time shifted signals on the solutions to these equations might be, I would rather not have to worry about it in the first place.

In matlab code this would be something like ifft(fft(ifft(fft(g)))) to get the original 'g' back.

you do get it back. but you have to perform the same shift you did after you got the bathtub.

in other words,

i got a function called x

X1 = fft(x) [indicies represent 0 to fs/2, -fs/2 to 0]
X2 = fftshift(fft(x)) [ indicies represent -fs/2 to fs/2]

after you do your processing with X2, you have to do another fftshift.
x = ifft(fftshift(X2))
ifft wants the [0 to fs/2, -fs to 0] version of your frequency response.

if you do ifft(fft(x) straight up you would get x in matlab.

also, when you said you put the gaussian at x = 3, you didn't actually center the fft of the gaussian at k = 3 did you? that would certainly screw things up.

 

[Biftheunderstudy]

Hmm, your right. If you do ifft(fft(x)) in matlab you get x back exactly If you do this with the c libraries you do not. Now to find out what Matlab does differently than the FFTW....

Edit:
Even stranger, the act of moving the fft routine calls into its own separate subroutine fixed the problem. Now I can do the equivalent of ifft(fft(x)) = x using the c libraries.

 

[TecHNooB]

Hmm, your right. If you do ifft(fft(x)) in matlab you get x back exactly If you do this with the c libraries you do not. Now to find out what Matlab does differently than the FFTW....

what exactly is different? is it shifted? or just scaled vertically?

 

[TecHNooB]

Edit:
Even stranger, the act of moving the fft routine calls into its own separate subroutine fixed the problem. Now I can do the equivalent of ifft(fft(x)) = x using the c libraries.

you should be well on your way now