Got Questions? Get Answers.
Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
fft on square wave

Subject: fft on square wave

From: Tomas

Date: 9 May, 2009 20:43:01

Message: 1 of 21

I want to compute the fourier coefficients of a square wave with height 1 and length 1 . According to the real fourier transform, the result should be F(z)=sin(pi*z)/(pi*z). This means the maximum value should be 1.

Can anyone explain me why I get growing max values of the discrete fourier coefficients as the smaller discretization step I use ?

Here is my code:

% "square function"
f = @(t) (heaviside(t+1/2) - heaviside(t-1/2));

% compute the fourier coefficients
t = linspace(-2,2);
F = fft(f(t));

% plot the coefficients (should be a real function)
plot(t,fftshift(real(F)))

 

Subject: fft on square wave

From: TideMan

Date: 9 May, 2009 21:12:21

Message: 2 of 21

On May 10, 8:43 am, "Tomas " <tomas...@gmail.com> wrote:
> I want to compute the fourier coefficients of a square wave with height 1=
 and length 1 . According to the real fourier transform, the result should =
be F(z)=sin(pi*z)/(pi*z). This means the maximum value should be 1.
>
> Can anyone explain me why I get growing max values of the discrete fourie=
r coefficients as the smaller discretization step I use ?
>
> Here is my code:
>
> % "square function"
> f = @(t) (heaviside(t+1/2) - heaviside(t-1/2));
>
> % compute the fourier coefficients
> t = linspace(-2,2);
> F = fft(f(t));
>
> % plot the coefficients (should be a real function)
> plot(t,fftshift(real(F)))

Can't you work it out?
When you halve the step, what happens to the magnitude of the DFT?
Now halve it again to check that the pattern is correct.

Hint: halving the step is equivalent to doubling the number of data

Subject: fft on square wave

From: Tomas

Date: 9 May, 2009 21:45:03

Message: 3 of 21

> Can't you work it out?
> When you halve the step, what happens to the magnitude of the DFT?
> Now halve it again to check that the pattern is correct.
>
> Hint: halving the step is equivalent to doubling the number of data

The point is that I don?t understand why the magnitude anyhow could be greater than 1 for the fourier transform of this square function. That is because the analytic answer is the sinc function (check wikipedia) and can possible not be greater than one.

Subject: fft on square wave

From: TideMan

Date: 9 May, 2009 22:08:15

Message: 4 of 21

On May 10, 9:45 am, "Tomas " <tomas...@gmail.com> wrote:
> > Can't you work it out?
> > When you halve the step, what happens to the magnitude of the DFT?
> > Now halve it again to check that the pattern is correct.
>
> > Hint: halving the step is equivalent to doubling the number of data
>
> The point is that I don?t understand why the magnitude anyhow could be gr=
eater than 1 for the fourier transform of this square function. That is bec=
ause the analytic answer is the sinc function (check wikipedia) and can pos=
sible not be greater than one.

Well, short of writing your own FFT, you'll just have to live with it
and get on with life.

Subject: fft on square wave

From: Matt

Date: 10 May, 2009 06:17:01

Message: 5 of 21

"Tomas " <tomasnic@gmail.com> wrote in message <gu4psl$oi$1@fred.mathworks.com>...
> I want to compute the fourier coefficients of a square wave with height 1 and length 1 . According to the real fourier transform, the result should be F(z)=sin(pi*z)/(pi*z). This means the maximum value should be 1.
>
> Can anyone explain me why I get growing max values of the discrete fourier coefficients as the smaller discretization step I use ?
-----

The fft does NOT approximate the Fourier Transform until you multiply it by your sampling interval. Remember, the continuous fourier transform is an integral while the DFT is a sum.

Incidentally, be aware that MATLAB's fft measures frequency in Hertz, not radians (see doc fft). Your expression for the sinc F(z)=sin(pi*z)/(pi*z) suggests you are looking for the latter.


> % compute the fourier coefficients
> t = linspace(-2,2);

This is not the way you want to generate your time axis. The fft() assumes that one of your time samples is t=0, which will not be the case with t=linspace(-2,2).

The formula for the axis you want to use is

( -ceil((N-1)/2):N-1-ceil((N-1)/2) )*SamplingInterval;
 %N is desired number of samples

With your current sampling, I'm pretty sure you're not getting a purely real spectrum.

Subject: fft on square wave

From: Tomas

Date: 10 May, 2009 10:24:01

Message: 6 of 21

> The fft does NOT approximate the Fourier Transform until you multiply it by your sampling interval. Remember, the continuous fourier transform is an integral while the DFT is a sum.

Can you explain this? According to the documentation the fft is a divide and conquer implementation of the discrete version of the real Fourier transform (integral replaced with a trapezoidal rule)

> Incidentally, be aware that MATLAB's fft measures frequency in Hertz, not radians (see doc fft). Your expression for the sinc F(z)=sin(pi*z)/(pi*z) suggests you are looking for the latter.

The measure depends on the function.

> This is not the way you want to generate your time axis. The fft() assumes that one of your time samples is t=0, which will not be the case with t=linspace(-2,2).

I used a new time discretization, containing t = 0:

% "square function"
f = @(t) (heaviside(t+1/2) - heaviside(t-1/2));

% compute the fourier coefficients
t = -2:0.01:2;
F = fft(f(t));

% plot the coefficients (should be a real function)
plot(t,fftshift(real(F)))

Now, the result ( F vector) is just non sense (containing NaN elements).
Can someone explain me this? Maybe im just really slow today (sunday morning).

Subject: fft on square wave

From: Matt

Date: 10 May, 2009 14:49:01

Message: 7 of 21

"Tomas " <tomasnic@gmail.com> wrote in message <gu6a01$5hq$1@fred.mathworks.com>...
> > The fft does NOT approximate the Fourier Transform until you multiply it by your sampling interval. Remember, the continuous fourier transform is an integral while the DFT is a sum.
>
> Can you explain this? According to the documentation the fft is a divide and conquer implementation of the discrete version of the real Fourier transform (integral replaced with a trapezoidal rule)
--------------------------------------

Don't know what documentation you're refering to. The fairly standard formula used by MATLAB's FFT (see "help fft") is

                     N
       X(k) = sum x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
                    n=1


It's just a sum of complex exponentials obtained by discretizing the integrand of the Fourier integral. It has no way of knowing how you are sampling time:


> I used a new time discretization, containing t = 0:
>
> % "square function"
> f = @(t) (heaviside(t+1/2) - heaviside(t-1/2));
>
> % compute the fourier coefficients
> t = -2:0.01:2;
> F = fft(f(t));

I think you really want

fft(f(ifftshift(t)))

because the fft() function expects that t(0)=0

> % plot the coefficients (should be a real function)
> plot(t,fftshift(real(F)))
>
> Now, the result ( F vector) is just non sense (containing NaN elements).
> Can someone explain me this? Maybe im just really slow today (sunday morning).

I can only suppose it's because of how you implement heaviside(). We can't see this because heaviside is not a native MATLAB function (or else it requires a ToolBox).


 

Subject: fft on square wave

From: Greg Heath

Date: 10 May, 2009 15:38:00

Message: 8 of 21

On May 9, 4:43 pm, "Tomas " <tomas...@gmail.com> wrote:
> I want to compute the fourier coefficients of a square wave with height 1=
 and length 1 . According to the real fourier transform, the result should =
be F(z)=sin(pi*z)/(pi*z). This means the maximum value should be 1.
>
> Can anyone explain me why I get growing max values of the discrete fourie=
r coefficients as the smaller discretization step I use ?
>
> Here is my code:
>
> % "square function"
> f = @(t) (heaviside(t+1/2) - heaviside(t-1/2));
>
> % compute the fourier coefficients
> t = linspace(-2,2);
> F = fft(f(t));
>
> % plot the coefficients (should be a real function)
> plot(t,fftshift(real(F)))

If by "real" you mean continuous time, note that
when you approximate an integal over t by a
uniformly sampled sum, there is a factor

dt = 1/Fs = T/N

that needs to be included.

Hope this helps.

Greg

Subject: fft on square wave

From: Tomas

Date: 10 May, 2009 18:30:03

Message: 9 of 21

I didn?t know the heaviside was not a native function. Here is an equal implementation of the square function:

function res = squarefunction(x)
    for i = 1 : length(x)
        if x(i) < -1 || x(i) > 1
            res(i) = 0;
        else
            res(i) = 1;
        end
    end
end

The Fourier transform of this function should in continous time be sin(t*pi)/(t*pi) and nowhere be greater than 1.

> If by "real" you mean continuous time, note that
> when you approximate an integal over t by a
> uniformly sampled sum, there is a factor
>
> dt = 1/Fs = T/N
>
> that needs to be included.

Yes, I ment "real" as in "real numbers". I shall from now use "continous time". Could you please explain the equation dt = 1/Fs = T/N, what the symbols Fs, T, N mean?

Thanks!

Subject: fft on square wave

From: Matt

Date: 11 May, 2009 00:25:03

Message: 10 of 21

"Tomas " <tomasnic@gmail.com> wrote in message <gu76fb$qv$1@fred.mathworks.com>...
> I didn?t know the heaviside was not a native function. Here is an equal implementation of the square function:
>
> function res = squarefunction(x)
> for i = 1 : length(x)
> if x(i) < -1 || x(i) > 1
> res(i) = 0;
> else
> res(i) = 1;
> end
> end
> end
>
> The Fourier transform of this function should in continous time be sin(t*pi)/(t*pi) and nowhere be greater than 1.
-------------------------

That won't help us identify the cause of the NaNs. It will only help if we see precisely what you are using.

Subject: fft on square wave

From: Matt

Date: 11 May, 2009 00:42:01

Message: 11 of 21

"Tomas " <tomasnic@gmail.com> wrote in message <gu6a01$5hq$1@fred.mathworks.com>...

> Now, the result ( F vector) is just non sense (containing NaN elements).
> Can someone explain me this? Maybe im just really slow today (sunday morning).


You should check whether there are NaNs in f(t),e.g. using

nnz(isnan(f(t)))

Subject: fft on square wave

From: Tomas

Date: 11 May, 2009 07:30:04

Message: 12 of 21

I found out the reason of the Nan. It?s because of the implementation of the heaviside function I used:

>> help heaviside
 HEAVISIDE Step function.
     HEAVISIDE(X) is 0 for X < 0, 1 for X > 0, and NaN for X == 0.
     HEAVISIDE(X) is not a function in the strict sense.
     See also dirac.

    Overloaded methods:
       sym/heaviside

    Reference page in Help browser
       doc heaviside

>>

Its returning NaN for X == 0.

I still hope that someone can explain me the equation Greg came with

dt = 1/Fs = T/N

I shall then multiply each fourier coefficient with T/N to prevent the "blow up" at t = 0 when I increase the number of time steps?

Subject: fft on square wave

From: Matt

Date: 11 May, 2009 12:48:01

Message: 13 of 21

"Tomas " <tomasnic@gmail.com> wrote in message <gu8k5s$knf$1@fred.mathworks.com>...

> I still hope that someone can explain me the equation Greg came with
>
> dt = 1/Fs = T/N
>
> I shall then multiply each fourier coefficient with T/N to prevent the "blow up" at t = 0 when I increase the number of time steps?

Greg is telling you the same thing I told you, which is to multiply the fft() output by your sampling interval, dt.

Beyond that, I don't quite follow his notation, but I think he once you to realize that once you've chosen a number of samples N, the following relationship holds between time sampling and frequency sampling:

N*TimeSamplingInterval*FrequencySamplingInterval=1

Subject: fft on square wave

From: Greg Heath

Date: 13 May, 2009 02:38:31

Message: 14 of 21

On May 9, 4:43 pm, "Tomas " <tomas...@gmail.com> wrote:
> I want to compute the fourier coefficients of a square wave
> with height 1 and length 1 . According to the real fourier
> transform, the result should be F(z)=sin(pi*z)/(pi*z). This
> means the maximum value should be 1.
>
> Can anyone explain me why I get growing max values of the
> discrete fourier coefficients as the smaller discretization
> step I use ?

You left out a factor dt = 1/Fs

> Here is my code:
>
> % "square function"
> f = @(t) (heaviside(t+1/2) - heaviside(t-1/2));

1. Version 6.5 does not have this function.
2. fft assumes sampling begins at t = 0

> % compute the fourier coefficients
> t = linspace(-2,2);

3. With 2 arguments, linspace always yields 100 points.

> F ÿt(f(t));
> % plot the coefficients (should be a real function)

Should never make that assumption always check
imag(F).

> plot(t,fftshift(real(F)))

There is nothing in this code that allows you to change
N (number of points).

Did you actually run this code?

%%%%%%%%%%%%%%%%%%%%%%%%%%

close all, clear all, k = 0
N = 16
for i = 1:6
    N = 2*N
    T = 2
    dt = T/N
    t = dt*(0:N-1);
    z = 0.5*(sign(t-1/2)-sign(t-3/2));
    k = k+1, figure(k)
    subplot(2,2,1)
    plot(t,z)
    axis([0 2 0 2])

    df = 1/T
    f = df*(0:N-1);
    Q = ceil((N+1)/2)
    fQ = (Q-1)*df
    fc = f-fQ;
    Z = fftshift(dt*fft(z));
    R = abs(Z);
    X = real(Z);
    Y = imag(Z);
    P = angle(Z);
    subplot(2,2,2)
    plot(fc,R)
    subplot(2,2,3)
    plot(fc,X)
    subplot(2,2,4)
    plot(fc,Y)

    n(i) = N;
    Xmax(i) = max(X);
    Ymax(i) = max(Y);

end
for j = k:-1:1
    figure(j)
end
format short g
summary = [ n' Xmax' Ymax']

Hope this helps.

Greg

Subject: fft on square wave

From: Greg Heath

Date: 13 May, 2009 03:08:03

Message: 15 of 21

On May 11, 3:30 am, "Tomas " <tomas...@gmail.com> wrote:
> I found out the reason of the Nan. It?s because of the implementation of =
the heaviside function I used:
>
> >> help heaviside
>
>  HEAVISIDE    Step function.
>      HEAVISIDE(X) is 0 for X < 0, 1 for X > 0, and NaN for X ===
 0.
>      HEAVISIDE(X) is not a function in the strict sense.
>      See also dirac.

I don't know where you got that from

Heaviside defined H(0) = H(0+) = lim(dx --> 0){ H(0+|dx|) } = 1
because he was interested in initial value problems using the
Laplace Transform whose inverse is defined on (0+,inf).
You will recall the initial value theorem of Laplace
Transform theory

lim(s-->0){ s*F(s) } = f(0+)

Since then, the doublesided Laplace Transform, Fourier
Transform and Fourier series have been in frequent
use. For these the value at a discontinuity is the average
of the left and right hand limits. Things turn out OK
if you define

H(t) = 0.5*(1+sign(t))

Then H(0) = 0.5,

>     Overloaded methods:
>        sym/heaviside
>
>     Reference page in Help browser
>        doc heaviside
>
> Its returning NaN for X == 0.

Someone needs to change that.

> I still hope that someone can explain me the equation Greg came with
>
> dt = 1/Fs = T/N
>
> I shall then multiply each fourier coefficient with T/N to prevent the "b=
low up" at t = 0 when I increase the number of time steps?

No, the "blow up at t=0 has to do with the bogus definition of the
heaviside function.

Fs = sampling frequency
dt = 1/Fs = time intervals between samples
T = N*dt = fft imposed period ( f(t+T) = f(t) )

See the comments and code in my recent post,

Hope this helps.

Greg

Subject: fft on square wave

From: Tomas

Date: 13 May, 2009 10:57:01

Message: 16 of 21

Thanks a lot for your explanation and the program, Greg! I now fully understand what went wrong.

Subject: fft on square wave

From: Matt

Date: 13 May, 2009 15:17:02

Message: 17 of 21

Greg Heath <heath@alumni.brown.edu> wrote in message <530b5762-c7ff-43c5-bd95-a10c347355bd@t11g2000vbc.googlegroups.com>...

> > F ?t(f(t));
> > % plot the coefficients (should be a real function)
>
> Should never make that assumption always check
> imag(F).
>

Tomas, since you're dealing with an even function (symmetric across t=0), the FFT will be real if you have an odd number N of samples and if t=0 is one of them.

Once you modified your time sampling axis to t=(-.2:.01:.2), this is now the case, so you can expect essentially real valued results. Of course, MATLAB doesn't know about the symmetry and will still return a complex result, but using imag(F) you should be able to verify that it contains zeroes up to numerical round-off.

Subject: fft on square wave

From: Dragon_Maja

Date: 24 Jun, 2009 12:30:15

Message: 18 of 21

Ciao I am new here.

I would like to ask for the hints of the problem:
I have a matrix Z=86x76, that I have got with confromal mapping like
zi=f(Z). And I have to find the time derivative (Z(t+dt)-Z(t))/dt in
computational coordinates. When I used the diff(Z) to find Z(t+dt)-Z
(t)) I have got the matrix 85x76, but it is not size that I need. Do
you have any hint how to find this time derivative.

Subject: fft on square wave

From: Dave Robinson

Date: 24 Jun, 2009 13:42:00

Message: 19 of 21

Dragon_Maja <marijavukicevic1@gmail.com> wrote in message <e786153c-8af8-4b18-8eb1-8fbecd574c35@q14g2000vbn.googlegroups.com>...
> Ciao I am new here.
>
> I would like to ask for the hints of the problem:
> I have a matrix Z=86x76, that I have got with confromal mapping like
> zi=f(Z). And I have to find the time derivative (Z(t+dt)-Z(t))/dt in
> computational coordinates. When I used the diff(Z) to find Z(t+dt)-Z
> (t)) I have got the matrix 85x76, but it is not size that I need. Do
> you have any hint how to find this time derivative.

You should have started a new thread - not added it to an existing one. Looking at it, it has nothing to do with "fft on square wave". What you might well find that guru's that can help you may well not even bother to read something on fft and squarewaves, and thus completely miss your post.

Whenever starting a new thread ALWAYS give it a meaningful title, so the people who are equipped to answer it will at least look at it, and not bypass it.

Regards

Dave Robinson

Subject: fft on square wave

From: Carlos Adrian Vargas Aguilera

Date: 24 Jun, 2009 21:44:01

Message: 20 of 21

WHAT A MESS!

You have three problems

1) frequency is never equal to time
2) MATLAB FFT theats the argument funtion as begining from zero. In your case the 0 is centered. Then, there is a time shift on f or a phase change on F to be considered, try:
3) you're a theoretician-man maybe MATHEMATICAtician (last but most important)

% --- COPY

% MODIFY:
T = 1; % square width
A = 1; % square amplitude
tmax = 2*T; % argument will go from -tmax to tmax
N = 100; % number of points

% "square function"
f = @(t) A*(heaviside(t+T/2) - heaviside(t-T/2));

% time from -tmax to tmax forced to pass through 0
DT = 2*tmax/100;
t = ((0:N-1)-floor(N/2))'*DT;

% fourier transform
w = t/(N*DT) / DT; % fix problem 2) (ugly way in fact)
F = fftshift(fft(f(t))*DT); % DT because of the "integration": dt

% moves phase because t(0)~=0 (fix problem 3), get it?)
F = F.*exp(1i*2*pi*w*t(1));

% sinc function
G = A*T*sinc(w*T); % this is your wiki

% plot the coefficients (should be a real function)
figure(1),clf
subplot(211), plot(t,f(t))
subplot(212), plot(w,G,'.-',w,[real(F) imag(F)])
legend('SINC','Re(FFT)','Im(FFT)')

% --- PASTE

So, who will fix the 3) problem... any?

Carlos Vargas

Subject: fft on square wave

From: Carlos Adrian Vargas Aguilera

Date: 24 Jun, 2009 22:45:03

Message: 21 of 21

DT = 2*tmax/100;
should be
DT = 2*tmax/N;

Well

Carlos Vargas

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us