Matlab fft and frequency vector

Hi,

I'm kind of newbie with matlab so pleace be patient. =) Well anyway, in one
algorithm I'm trying to implement with Matlab I need to divide a Fourier
transform with i*w (i is imaginary unit and w omega or 2*pi*f). So it's
something like this:

Fs = 501;   %sampling frequency
t=[-1:1/Fs:1];
x = cos(2*pi*10*t);

Xf = fft(x);

f = [-0.5*Fs:Fs/N:0.5*Fs-1/Fs]; frequency vector

xt1 = ifft(ifftshift(fftshift(Xf)./(j*2*pi*f)));

subplot(2,1,1);
plot(t,x);
grid on;

subplot(2,1,2);
plot(t,xt1);
grid on;

This works a little bit but at least the scaling goes wrong. And I've also
noticed there can't be any zeros in the frequency vector, because then
everithing goes to infinity. I know that diveding with i*w in the
frequency-domain is same as integration in time-domain, but I don't really
know how to do integration with Matlab. So if you have any correction ideas
or some totally different kind of approach in mind, help would be greatly
appreciated.

-Mika

"Mika" <ee8460@kolumbus.fi> wrote in message 
news:dvuq41$1ll$1@phys-news4.kolumbus.fi...
> Hi,
>
> I'm kind of newbie with matlab so pleace be patient. =) Well anyway, in 
> one
> algorithm I'm trying to implement with Matlab I need to divide a Fourier
> transform with i*w (i is imaginary unit and w omega or 2*pi*f). So it's
> something like this:
>
> Fs = 501;   %sampling frequency
> t=[-1:1/Fs:1];
> x = cos(2*pi*10*t);
>
> Xf = fft(x);
>
> f = [-0.5*Fs:Fs/N:0.5*Fs-1/Fs]; frequency vector
>
> xt1 = ifft(ifftshift(fftshift(Xf)./(j*2*pi*f)));
>
> subplot(2,1,1);
> plot(t,x);
> grid on;
>
> subplot(2,1,2);
> plot(t,xt1);
> grid on;
>
> This works a little bit but at least the scaling goes wrong. And I've also
> noticed there can't be any zeros in the frequency vector, because then
> everithing goes to infinity. I know that diveding with i*w in the
> frequency-domain is same as integration in time-domain, but I don't really
> know how to do integration with Matlab. So if you have any correction 
> ideas
> or some totally different kind of approach in mind, help would be greatly
> appreciated.
>
> -Mika
>
> 

"Mika" <ee8460@kolumbus.fi> wrote in message 
news:dvuq41$1ll$1@phys-news4.kolumbus.fi...
> Hi,
>
> I'm kind of newbie with matlab so pleace be patient. =) Well anyway, in 
> one
> algorithm I'm trying to implement with Matlab I need to divide a Fourier
> transform with i*w (i is imaginary unit and w omega or 2*pi*f). So it's
> something like this:
>
> Fs = 501;   %sampling frequency
> t=[-1:1/Fs:1];
> x = cos(2*pi*10*t);
>
> Xf = fft(x);
>
> f = [-0.5*Fs:Fs/N:0.5*Fs-1/Fs]; frequency vector
>
> xt1 = ifft(ifftshift(fftshift(Xf)./(j*2*pi*f)));
>
> subplot(2,1,1);
> plot(t,x);
> grid on;
>
> subplot(2,1,2);
> plot(t,xt1);
> grid on;
>
> This works a little bit but at least the scaling goes wrong. And I've also
> noticed there can't be any zeros in the frequency vector, because then
> everithing goes to infinity. I know that diveding with i*w in the
> frequency-domain is same as integration in time-domain, but I don't really
> know how to do integration with Matlab. So if you have any correction 
> ideas
> or some totally different kind of approach in mind, help would be greatly
> appreciated.
>
> -Mika

I don't know, it looks like you're doing quite a few unecessary fft/ifft 
operations!

Xf > fft of x with normal registration
fftshift(Xf) hmmmmm... I have NO idea what you're attempting here.
  what you will get is a sort of ifft of Xf that is shifted around in time.
ifftshift(the above) what you get is a sort of fft of x that is shifted 
around in frequency
ifft (of the above) what you get is an ifft of that strange fft of that 
strange ifft of that fft
Whew! that's enough to make my mind spin!
I have no idea what that would be.

Are you sure you don't want something like this:

fftshift(x)=Xf
then weigh by the intended weighting function Xf > Xfw
ifftshift (Xfw)
done!
I'm not sure but I presume that fftshift followed by ifftshift gets you back 
to x.

Fred 

> I don't know, it looks like you're doing quite a few unecessary fft/ifft
> operations!
>

Yeah you might be right =)

> Xf > fft of x with normal registration
> fftshift(Xf) hmmmmm... I have NO idea what you're attempting here.
>   what you will get is a sort of ifft of Xf that is shifted around in
time.
> ifftshift(the above) what you get is a sort of fft of x that is shifted
> around in frequency
> ifft (of the above) what you get is an ifft of that strange fft of that
> strange ifft of that fft
> Whew! that's enough to make my mind spin!
> I have no idea what that would be.
>

fftshift means that it shifts zero-frequency component to center of
spectrum. Because Matlab's fft puts the zero frequency to the beginnig. And
ifftshift undoes the effects of fftshift. But they are not really necessary
because you can get the same results using different frequency vector. Like
this

f1 = [1/Fs:Fs/N:0.5*Fs];
f2 = [-0.5*Fs:Fs/N:0-Fs/N];
f = [f1 f2];
xt1 = ifft(fft(Xf)./(j*2*pi*f));

is equivalent with

f = [-0.5*Fs:Fs/N:0.5*Fs-1/Fs];
xt1 = ifft(ifftshift(fftshift(fft(Xf))./(j*2*pi*f)));

in the first one the frequency vector goes from 0 to Fs/2 and then
from -Fs/2 to 0. And in the second one the frequency vector goes from -Fs/2
to Fs/2. And my real problem is how I should choose the right frequency
vector?

> Are you sure you don't want something like this:
>
> fftshift(x)=Xf
> then weigh by the intended weighting function Xf > Xfw
> ifftshift (Xfw)
> done!
> I'm not sure but I presume that fftshift followed by ifftshift gets you
back
> to x.

Actually it should go something like
fft(x) = Xf
then the weighting Xf > Xfw
ifft(Xfw)

I just don't know the right weighting function. =)

-Mika

"Mika" <ee8460@kolumbus.fi> wrote in message 
news:dvv04a$50b$1@phys-news4.kolumbus.fi...


>
> I just don't know the right weighting function. =)
>
> -Mika

OK - what are you trying to accomplish?

Fred 

"Fred Marshall" <fmarshallx@remove_the_x.acm.org> wrote in message
news:jtWdncwv3qwoy77ZRVn-jQ@centurytel.net...
>
> "Mika" <ee8460@kolumbus.fi> wrote in message
> news:dvv04a$50b$1@phys-news4.kolumbus.fi...
>
>
> >
> > I just don't know the right weighting function. =)
> >
> > -Mika
>
> OK - what are you trying to accomplish?
>
I'm trying to count Lateral Energy Fraction (LEF) of room with the method
suggested by Kleiner 1989. In this method you measure the room impulse
response with two omni michrophones and from those two impulse responses
according to Kleiner you should get figure8 microphone response like this:

fig8 = ifft(fft(omni1-omni2)*c/(j*2*pi*f*b));

where omni1 and omni2 are the measured impulse responses, c speed of sound
and b distance between the microphones. But my problem is how I create this
frequency vector f in Matlab. I have tried things like

f = [-0.5*Fs:Fs/N:0.5*Fs-1/Fs]; and

f1 = [1/Fs:Fs/N:0.5*Fs];
f2 = [-0.5*Fs:Fs/N:0-Fs/N];
f = [f1 f2];  and

f = [0:Fs/N:Fs-Fs/N];

where Fs is the samplig frequency and N length of the impulse response.
However none of these seem to work.

-Mika

"Mika" <ee8460@kolumbus.fi> wrote in message 
news:e00fcj$bjm$1@phys-news4.kolumbus.fi...
>
> "Fred Marshall" <fmarshallx@remove_the_x.acm.org> wrote in message
> news:jtWdncwv3qwoy77ZRVn-jQ@centurytel.net...
>>
>> "Mika" <ee8460@kolumbus.fi> wrote in message
>> news:dvv04a$50b$1@phys-news4.kolumbus.fi...
>>
>>
>> >
>> > I just don't know the right weighting function. =)
>> >
>> > -Mika
>>
>> OK - what are you trying to accomplish?
>>
> I'm trying to count Lateral Energy Fraction (LEF) of room with the method
> suggested by Kleiner 1989. In this method you measure the room impulse
> response with two omni michrophones and from those two impulse responses
> according to Kleiner you should get figure8 microphone response like this:
>
> fig8 = ifft(fft(omni1-omni2)*c/(j*2*pi*f*b));
>
> where omni1 and omni2 are the measured impulse responses, c speed of sound
> and b distance between the microphones. But my problem is how I create 
> this
> frequency vector f in Matlab. I have tried things like
>
> f = [-0.5*Fs:Fs/N:0.5*Fs-1/Fs]; and
>
> f1 = [1/Fs:Fs/N:0.5*Fs];
> f2 = [-0.5*Fs:Fs/N:0-Fs/N];
> f = [f1 f2];  and
>
> f = [0:Fs/N:Fs-Fs/N];
>
> where Fs is the samplig frequency and N length of the impulse response.
> However none of these seem to work.
>
> -Mika

I don't know about LEF.
Your words suggest that a "microphone response" that is a beam pattern but 
I'm not sure.
"figure 8" surely suggests that...
So does the distance between microphones and the speed of sound being in 
there.

Accordingly it's unclear to me where the "fs" comes in to any of this.  For 
a beam pattern I'd probably use something like -pi to pi.

Now, if you're doing numeric calculations then indeed things need to be 
sampled or discrete.  So maybe the solution to your problem is to decide 
what the independent variables are and go from there.

In the response fig8, I see "f" as frequency.
If it's for a beam pattern then there's a pattern for *each* frequency but I 
don't see the spatial angles like theta or phi in the expression.  Maybe 
it's a single measure of each pattern?

Anyway, I'd be likely to compute the measure as a function of frequency 
going from zero to however high you need to go and not include negative 
frequencies as this appears to not be a signal but, rather, a measure.

But I'm guessing.

Fred 

"Mika" <ee8460@kolumbus.fi> writes:

> Hi,
>
> I'm kind of newbie with matlab so pleace be patient. =) Well anyway, in one
> algorithm I'm trying to implement with Matlab I need to divide a Fourier
> transform with i*w (i is imaginary unit and w omega or 2*pi*f). So it's
> something like this:
>
> Fs = 501;   %sampling frequency
> t=[-1:1/Fs:1];
> x = cos(2*pi*10*t);
>
> Xf = fft(x);
>
> f = [-0.5*Fs:Fs/N:0.5*Fs-1/Fs]; frequency vector
>
> xt1 = ifft(ifftshift(fftshift(Xf)./(j*2*pi*f)));
>
> subplot(2,1,1);
> plot(t,x);
> grid on;
>
> subplot(2,1,2);
> plot(t,xt1);
> grid on;
>
> This works a little bit but at least the scaling goes wrong. And I've also
> noticed there can't be any zeros in the frequency vector, because then
> everithing goes to infinity. I know that diveding with i*w in the
> frequency-domain is same as integration in time-domain, but I don't really
> know how to do integration with Matlab. So if you have any correction ideas
> or some totally different kind of approach in mind, help would be greatly
> appreciated.
>
> -Mika

Hi Mika,

Several points:

  1. Using the N-point DFT implies that the input signal x[n] is periodic
  with period N. Is this the assumption you wish to make?

  2. The integration property you describe has the following form for the
  discrete-time Fourier series (i.e., the DFT):

    \sum_{k=-\infty}^{n}  <==>  (1 / (1 - exp(-jk(2*pi/N)))) * a_k,

  where a_k is the kth DFT coefficient of x[n], 0 <= k <= N - 1.

  3. The fftshift/ifftshift operations are incorrect according to the 
  above property. 

  4. The division is not by j*2*pi*f, as for the continuous-time
  Fourier transform, but by 1 - exp(-jk(2*pi/N)), per the property
  above.

  5. The result is finite-valued only if there is no DC term in
  the Fourier series, i.e., a_0 = 0.

This is from "Signals and Systems," Oppenheim et al.

So, your Matlab code is fine - your theory is wrong. In general, you
cannot take a continuous-time Fourier transform property, discretize
it, and apply it to the DFT. 
-- 
%  Randy Yates                  % "Remember the good old 1980's, when 
%% Fuquay-Varina, NC            %  things were so uncomplicated?"
%%% 919-577-9882                % 'Ticket To The Moon' 
%%%% <yates@ieee.org>           % *Time*, Electric Light Orchestra
http://home.earthlink.net/~yatescr