IFFT problem rephrased...

I want to simulate the Inverse Fourier Transform of a brickwall/rectangle
function from -1500Hz to 1500Hz out of [-5000Hz, 5000Hz], the height is
1/Fd, where Fd=3000.

Theoratically I know the IFT should be sinc(Fd*t).

----------------------------------------------------
%Inverse FFT ...
Fs=10;           %this is the sampling rate for inverse FFT.
Fn=5000;         %Nyquist Rate.
Fd=3000;         %from -1500Hz to 1500Hz.
rect=inline('double(abs(f)<=0.5)', 'f');
f=[-Fn:1/Fs:Fn]; %frequency range from -5000Hz to +5000Hz
X=1/Fd*rect(f/Fd); %rect function between [-1500Hz, +1500Hz]&#4294967295;&#4294967295; height 1/Fd.

%sinc(Fd*t) <-> 1/Fd*rect(f/Fd); %theoratically

L=2^(round(log2(length(f))));

NewFs=L/Fs;
stepsize=1/NewFs;

%I don't know what should be the scaling
%factor in this ifft process...here I put Fs, but it was wrong...
x=ifft(ifftshift(X), L)*Fs;
plot([0:L-1]*Fs/L-Fs/2, abs(fftshift(x)));

------------------------------------------------------

The problems are:

1. I don't know what should be the scaling factor in the IFFT expression...
2. I don't know what range in time domain does the output vector x[1..L]
represent?
3. If I want to map the output range in time domain to any arbitrary range,
say [-2s to 2s], can I do that by adjusting the sampling rate in frequency
domain and the length of block L?
4. I really am confused by the fact that I need to put an additional
"fftshift()" on the IFFT output vector in order to reorder... it does not
make sense. There should have no such need to reorder...

Thank you very much!

"lucy" <losemind@yahoo.com> wrote in
news:cfuvad$r30$1@news.Stanford.EDU: 

> I want to simulate the Inverse Fourier Transform of a
> brickwall/rectangle function from -1500Hz to 1500Hz out of [-5000Hz,
> 5000Hz], the height is 1/Fd, where Fd=3000.
> 
> Theoratically I know the IFT should be sinc(Fd*t).
> 
> ----------------------------------------------------
> %Inverse FFT ...
> Fs=10;           %this is the sampling rate for inverse FFT.
> Fn=5000;         %Nyquist Rate.

The Nyquist rate is half the sampling rate, by definition. There is a 
real problem with something here

> Fd=3000;         %from -1500Hz to 1500Hz.
> rect=inline('double(abs(f)<=0.5)', 'f');
> f=[-Fn:1/Fs:Fn]; %frequency range from -5000Hz to +5000Hz
> X=1/Fd*rect(f/Fd); %rect function between [-1500Hz, +1500Hz]&#4294967295;&#4294967295; height
> 1/Fd. 
> 
> %sinc(Fd*t) <-> 1/Fd*rect(f/Fd); %theoratically
> 
> L=2^(round(log2(length(f))));
> 
> NewFs=L/Fs;
> stepsize=1/NewFs;
> 
> %I don't know what should be the scaling
> %factor in this ifft process...here I put Fs, but it was wrong...
> x=ifft(ifftshift(X), L)*Fs;
> plot([0:L-1]*Fs/L-Fs/2, abs(fftshift(x)));
> 
> ------------------------------------------------------
> 
> The problems are:
> 
> 1. I don't know what should be the scaling factor in the IFFT
> expression... 2. I don't know what range in time domain does the
> output vector x[1..L] represent?
> 3. If I want to map the output range in time domain to any arbitrary
> range, say [-2s to 2s], can I do that by adjusting the sampling rate
> in frequency domain and the length of block L?
> 4. I really am confused by the fact that I need to put an additional
> "fftshift()" on the IFFT output vector in order to reorder... it does
> not make sense. There should have no such need to reorder...
> 
> Thank you very much!
> 
> 
> 
> 
If "width" is the number of points in the sinc function you're trying 
for, and "sampfreq" is your sampling rate (which defines the range of 
your FFT), and "cutoff" is your desired cutoff frequency (the end of your 
rectangle) then 

ntemp=-1*fix(width/2):fix(width/2);
% compute impulse response of ideal lpf
% (that's the sinc function you're looking for!!)
impresp=cutoff*2/sampfrq*sinc(cutoff*2*pi*ntemp/pi/sampfrq);

Note that this is used here to design a filter of variable width.  I 
think in your case the width simply needs to be NFFT.

I'd really recommend picking up a copy of Oppenhiem and Schafer if you 
find yourself doing DSP.  You'd save yourself alot of time, and you'd 
understand the area better in the end.

There are some issues with Matlab's definition of the FFT not being 
scaled by a factor of 1/N, and that tends to bubble through in certain 
places, thus matlab results might differ from the textbook equations.  
This has been discussed ad nauseum in this newsgroup.  Compare the matlab 
definitions available in the docs to your textbook equations.  
Empirically, look at your results and figure out if you need to multiply 
by N or 1/N to make things work out.

Scott

Hi Lucy 
	What you should get is a sampled sinc (actually its a sampled dirichlet
function but why quibble). If you were to create a rectangle that
spanned the entire frequency range and performed your IDFT you would get
one sample that lands on the main lobe of the sinc and the rest land on
zero. That's not going to look like a sinc function to your eyeballs. If
you made your rectangle half of the frequency domain  span (and N twice
as big) you would get the same thing with the addition of one more
sample between each of the original samples. That might look more like
what you might expect. But if I understand correctly, you  are making
your rectangle at 30 percent of the frequency span so your sampled sinc
is probably not going to be easily recognized as such.

-jim 

lucy wrote:
> 
> I want to simulate the Inverse Fourier Transform of a brickwall/rectangle
> function from -1500Hz to 1500Hz out of [-5000Hz, 5000Hz], the height is
> 1/Fd, where Fd=3000.
> 
> Theoratically I know the IFT should be sinc(Fd*t).
> 
> ----------------------------------------------------
> %Inverse FFT ...
> Fs=10;           %this is the sampling rate for inverse FFT.
> Fn=5000;         %Nyquist Rate.
> Fd=3000;         %from -1500Hz to 1500Hz.
> rect=inline('double(abs(f)<=0.5)', 'f');
> f=[-Fn:1/Fs:Fn]; %frequency range from -5000Hz to +5000Hz
> X=1/Fd*rect(f/Fd); %rect function between [-1500Hz, +1500Hz]&#4294967295;&#4294967295; height 1/Fd.
> 
> %sinc(Fd*t) <-> 1/Fd*rect(f/Fd); %theoratically
> 
> L=2^(round(log2(length(f))));
> 
> NewFs=L/Fs;
> stepsize=1/NewFs;
> 
> %I don't know what should be the scaling
> %factor in this ifft process...here I put Fs, but it was wrong...
> x=ifft(ifftshift(X), L)*Fs;
> plot([0:L-1]*Fs/L-Fs/2, abs(fftshift(x)));
> 
> ------------------------------------------------------
> 
> The problems are:
> 
> 1. I don't know what should be the scaling factor in the IFFT expression...
> 2. I don't know what range in time domain does the output vector x[1..L]
> represent?
> 3. If I want to map the output range in time domain to any arbitrary range,
> say [-2s to 2s], can I do that by adjusting the sampling rate in frequency
> domain and the length of block L?
> 4. I really am confused by the fact that I need to put an additional
> "fftshift()" on the IFFT output vector in order to reorder... it does not
> make sense. There should have no such need to reorder...
> 
> Thank you very much!


-----= Posted via Newsfeeds.Com, Uncensored Usenet News =-----
http://www.newsfeeds.com - The #1 Newsgroup Service in the World!
-----==  Over 100,000 Newsgroups - 19 Different Servers! =-----

"Scott Seidman" <namdiesttocs@mindspring.com> wrote in message
news:Xns954954CF2B359scottseidmanmindspri@130.133.1.4...
> "lucy" <losemind@yahoo.com> wrote in
> news:cfuvad$r30$1@news.Stanford.EDU:
>
> > I want to simulate the Inverse Fourier Transform of a
> > brickwall/rectangle function from -1500Hz to 1500Hz out of [-5000Hz,
> > 5000Hz], the height is 1/Fd, where Fd=3000.
> >
> > Theoratically I know the IFT should be sinc(Fd*t).
> >
> > ----------------------------------------------------
> > %Inverse FFT ...
> > Fs=10;           %this is the sampling rate for inverse FFT.
> > Fn=5000;         %Nyquist Rate.
>
> The Nyquist rate is half the sampling rate, by definition. There is a
> real problem with something here
>
> > Fd=3000;         %from -1500Hz to 1500Hz.
> > rect=inline('double(abs(f)<=0.5)', 'f');
> > f=[-Fn:1/Fs:Fn]; %frequency range from -5000Hz to +5000Hz
> > X=1/Fd*rect(f/Fd); %rect function between [-1500Hz, +1500Hz]&#4294967295;&#4294967295; height
> > 1/Fd.
> >
> > %sinc(Fd*t) <-> 1/Fd*rect(f/Fd); %theoratically
> >
> > L=2^(round(log2(length(f))));
> >
> > NewFs=L/Fs;
> > stepsize=1/NewFs;
> >
> > %I don't know what should be the scaling
> > %factor in this ifft process...here I put Fs, but it was wrong...
> > x=ifft(ifftshift(X), L)*Fs;
> > plot([0:L-1]*Fs/L-Fs/2, abs(fftshift(x)));
> >
> > ------------------------------------------------------
> >
> > The problems are:
> >
> > 1. I don't know what should be the scaling factor in the IFFT
> > expression... 2. I don't know what range in time domain does the
> > output vector x[1..L] represent?
> > 3. If I want to map the output range in time domain to any arbitrary
> > range, say [-2s to 2s], can I do that by adjusting the sampling rate
> > in frequency domain and the length of block L?
> > 4. I really am confused by the fact that I need to put an additional
> > "fftshift()" on the IFFT output vector in order to reorder... it does
> > not make sense. There should have no such need to reorder...
> >
> > Thank you very much!
> >
> >
> >
> >
> If "width" is the number of points in the sinc function you're trying
> for, and "sampfreq" is your sampling rate (which defines the range of
> your FFT), and "cutoff" is your desired cutoff frequency (the end of your
> rectangle) then
>
> ntemp=-1*fix(width/2):fix(width/2);
> % compute impulse response of ideal lpf
> % (that's the sinc function you're looking for!!)
> impresp=cutoff*2/sampfrq*sinc(cutoff*2*pi*ntemp/pi/sampfrq);
>
> Note that this is used here to design a filter of variable width.  I
> think in your case the width simply needs to be NFFT.
>


Hi Scott,

Thank you very much for your help. I have read those O & S DSP books on a
self-study basis... now I am trying this kind of simple stuff. The books
read all OK! But when you turn from books to Matlab, I found there are so
much strange things that don't make sense...

For example:

plot([0:L-1]*Fs/L-Fs/2, abs(fftshift(x)));

if I don't put "abs" here, the plot is wierd, but if I put abs here,

all negative values are gone...

But sinc() function should have negative values...

What do you think?

"lucy" <losemind@yahoo.com> wrote in news:cg054o$k2b$1@news.Stanford.EDU:

> Hi Scott,
> 
> Thank you very much for your help. I have read those O & S DSP books on a
> self-study basis... now I am trying this kind of simple stuff. The books
> read all OK! But when you turn from books to Matlab, I found there are so
> much strange things that don't make sense...
> 
> For example:
> 
> plot([0:L-1]*Fs/L-Fs/2, abs(fftshift(x)));
> 
> if I don't put "abs" here, the plot is wierd, but if I put abs here,
> 
> all negative values are gone...
> 
> But sinc() function should have negative values...
> 
> What do you think?
> 
> 
> 

fft appropriately produces complex output, with a real and an imaginary 
part.  ifft tends to do the same sometimes, but usually because of roundoff 
errors.  For the ifft case, take the REAL part of the answer. For the fft 
case, the whole story is told by the abs(fftshift(...)) and angle(fftshift
(...)).  If you don't follow how the magnitude and angle describe the fft 
output, then you need to remediate...its beyond what you can expect out of 
a forum like this.

Scott

jim <"N0sp"@m.sjedging@mwt.net> wrote...
> 	What you should get is a sampled sinc (actually its a sampled dirichlet
> function but why quibble).

The *discrete* Fourier transform of a rectangular window function is
neither a sampled sinc nor a sampled Dirichlet function (if you're
thinking of the same Dirichlet function I know of, it only takes on
two values).

For example, suppose you have an N-point DFT of a function that is 1
for 0..n-1 and 0 elsewhere, the DFT (or IDFT, modulo complex
conjgation) is (you can simply sum the geometric series):

          exp(-pi*i*(n-1)*k/N) * sin(pi*n*k/2) / sin(pi*k/2)

The exp(...) is just a phase shift from the fact that the window is
not centered at 0.  The remainder is *not* a sinc, but is a sin/sin
function (which goes to a sinc function as N -> infinity (for fixed
frequency and window width).

Cordially,
Steven G. Johnson

"Scott Seidman" <namdiesttocs@mindspring.com> wrote in message
news:Xns954992AB49F80scottseidmanmindspri@130.133.1.4...
> "lucy" <losemind@yahoo.com> wrote in news:cg054o$k2b$1@news.Stanford.EDU:
>
> > Hi Scott,
> >
> > Thank you very much for your help. I have read those O & S DSP books on
a
> > self-study basis... now I am trying this kind of simple stuff. The books
> > read all OK! But when you turn from books to Matlab, I found there are
so
> > much strange things that don't make sense...
> >
> > For example:
> >
> > plot([0:L-1]*Fs/L-Fs/2, abs(fftshift(x)));
> >
> > if I don't put "abs" here, the plot is wierd, but if I put abs here,
> >
> > all negative values are gone...
> >
> > But sinc() function should have negative values...
> >
> > What do you think?
> >
> >
> >
>
> fft appropriately produces complex output, with a real and an imaginary
> part.  ifft tends to do the same sometimes, but usually because of
roundoff
> errors.  For the ifft case, take the REAL part of the answer. For the fft
> case, the whole story is told by the abs(fftshift(...)) and angle(fftshift
> (...)).  If you don't follow how the magnitude and angle describe the fft
> output, then you need to remediate...its beyond what you can expect out of
> a forum like this.
>
> Scott


Hi Scott,

Thank you very much. But I was talking about the IFFT...

Here is the snippet of the code in my original post.

x=ifft(ifftshift(X), L)*Fs;
plot([0:L-1]*Fs/L-Fs/2, abs(fftshift(x)));

If I don't use "abs" on the result of "ifft", the result is very wierd,

but if I apply "abs" on the result of "ifft"(as shown above),

I lost all negative components...

What's wrong with this?

"lucy" <losemind@yahoo.com> wrote in message
news:cg09pd$n5d$1@news.Stanford.EDU...
>
> "Scott Seidman" <namdiesttocs@mindspring.com> wrote in message
> news:Xns954992AB49F80scottseidmanmindspri@130.133.1.4...
> > "lucy" <losemind@yahoo.com> wrote in
news:cg054o$k2b$1@news.Stanford.EDU:
> >
> > > Hi Scott,
> > >
> > > Thank you very much for your help. I have read those O & S DSP books
on
> a
> > > self-study basis... now I am trying this kind of simple stuff. The
books
> > > read all OK! But when you turn from books to Matlab, I found there are
> so
> > > much strange things that don't make sense...
> > >
> > > For example:
> > >
> > > plot([0:L-1]*Fs/L-Fs/2, abs(fftshift(x)));
> > >
> > > if I don't put "abs" here, the plot is wierd, but if I put abs here,
> > >
> > > all negative values are gone...
> > >
> > > But sinc() function should have negative values...
> > >
> > > What do you think?
> > >
> > >
> > >
> >
> > fft appropriately produces complex output, with a real and an imaginary
> > part.  ifft tends to do the same sometimes, but usually because of
> roundoff
> > errors.  For the ifft case, take the REAL part of the answer. For the
fft
> > case, the whole story is told by the abs(fftshift(...)) and
angle(fftshift
> > (...)).  If you don't follow how the magnitude and angle describe the
fft
> > output, then you need to remediate...its beyond what you can expect out
of
> > a forum like this.
> >
> > Scott
>
>
> Hi Scott,
>
> Thank you very much. But I was talking about the IFFT...
>
> Here is the snippet of the code in my original post.
>
> x=ifft(ifftshift(X), L)*Fs;
> plot([0:L-1]*Fs/L-Fs/2, abs(fftshift(x)));
>
> If I don't use "abs" on the result of "ifft", the result is very wierd,

That's because the output of the ifft is a complex number (with a real part
and an imaginary part) and you are trying to plot the real/imaginary pairs.

> but if I apply "abs" on the result of "ifft"(as shown above),
>
> I lost all negative components...

Try 'real(fftshift(x))' instead.

Cheers
Bhaskar

Steven G. Johnson wrote:
> jim <"N0sp"@m.sjedging@mwt.net> wrote...
> 
>>	What you should get is a sampled sinc (actually its a sampled dirichlet
>>function but why quibble).
> 
> 
> The *discrete* Fourier transform of a rectangular window function is
> neither a sampled sinc nor a sampled Dirichlet function (if you're
> thinking of the same Dirichlet function I know of, it only takes on
> two values).
> 
> For example, suppose you have an N-point DFT of a function that is 1
> for 0..n-1 and 0 elsewhere, the DFT (or IDFT, modulo complex
> conjgation) is (you can simply sum the geometric series):
> 
>           exp(-pi*i*(n-1)*k/N) * sin(pi*n*k/2) / sin(pi*k/2)
> 
> The exp(...) is just a phase shift from the fact that the window is
> not centered at 0.  The remainder is *not* a sinc, but is a sin/sin
> function (which goes to a sinc function as N -> infinity (for fixed
> frequency and window width).
> 
> Cordially,
> Steven G. Johnson

In Matlab, the function diric from the signal processing toolbox is 
called the Dirichlet function and given as:

   sin(N x /2)/[N sim(x/2)]

I've used it to calculate beam patterns for uniform linear arrays.


What's your Dirichlet function look like?

Stan Pawlukiewicz wrote:
> 
> Steven G. Johnson wrote:
> > jim <"N0sp"@m.sjedging@mwt.net> wrote...
> >
> >>      What you should get is a sampled sinc (actually its a sampled dirichlet
> >>function but why quibble).
> >
> >
> > The *discrete* Fourier transform of a rectangular window function is
> > neither a sampled sinc nor a sampled Dirichlet function (if you're
> > thinking of the same Dirichlet function I know of, it only takes on
> > two values).
> >
> > For example, suppose you have an N-point DFT of a function that is 1
> > for 0..n-1 and 0 elsewhere, the DFT (or IDFT, modulo complex
> > conjgation) is (you can simply sum the geometric series):
> >
> >           exp(-pi*i*(n-1)*k/N) * sin(pi*n*k/2) / sin(pi*k/2)
> >
> > The exp(...) is just a phase shift from the fact that the window is
> > not centered at 0.  The remainder is *not* a sinc, but is a sin/sin
> > function (which goes to a sinc function as N -> infinity (for fixed
> > frequency and window width).
> >
> > Cordially,
> > Steven G. Johnson
> 
> In Matlab, the function diric from the signal processing toolbox is
> called the Dirichlet function and given as:
> 
>    sin(N x /2)/[N sin(x/2)]
> 
> I've used it to calculate beam patterns for uniform linear arrays.
> 
> What's your Dirichlet function look like?

My newsreader is not showing Steven's reply so thank-you for replying so
that I might see it.

	I assume Steven is saying that it is shifted and scaled in x so that
makes it not a Dirichlet function. I don't know, sounds like a quibble
to me. The dirichlet function can be derived as he said by summing the
geometric series so what's the difference.

-jim


-----= Posted via Newsfeeds.Com, Uncensored Usenet News =-----
http://www.newsfeeds.com - The #1 Newsgroup Service in the World!
-----==  Over 100,000 Newsgroups - 19 Different Servers! =-----