CODEWITHSUNDEEP
fftmatlab
Richard
Richard

Reputation: 3296

How to compute Fourier coefficients with MATLAB

I'm trying to compute the Fourier coefficients for a waveform using MATLAB. The coefficients can be computed using the following formulas:

enter image description here

enter image description here

T is chosen to be 1 which gives omega = 2pi.

However I'm having issues performing the integrals. The functions are are triangle wave (Which can be generated using sawtooth(t,0.5) if I'm not mistaking) as well as a square wave.

I've tried with the following code (For the triangle wave):

    function [ a0,am,bm ] = test( numTerms )
            b_m = zeros(1,numTerms);
            w=2*pi;
            for i = 1:numTerms
                    f1 = @(t) sawtooth(t,0.5).*cos(i*w*t);
                    f2 = @(t) sawtooth(t,0.5).*sin(i*w*t);
                    am(i) = 2*quad(f1,0,1);
                    bm(i) = 2*quad(f2,0,1);
            end
    end

However it's not getting anywhere near the values I need. The b_m coefficients are given for a triangle wave and are supposed to be 1/m^2 and -1/m^2 when m is odd alternating beginning with the positive term.

The major issue for me is that I don't quite understand how integrals work in MATLAB and I'm not sure whether or not the approach I've chosen works.

Edit: To clairify, this is the form that I'm looking to write the function on when the coefficients have been determined:

enter image description here

Here's an attempt using fft:

    function [ a0,am,bm ] = test( numTerms )
            T=2*pi;
            w=1;
            t = [0:0.1:2];
            f = fft(sawtooth(t,0.5));
            am = real(f);
            bm = imag(f);
            func = num2str(f(1));
            for i = 1:numTerms
                    func = strcat(func,'+',num2str(am(i)),'*cos(',num2str(i*w),'*t)','+',num2str(bm(i)),'*sin(',num2str(i*w),'*t)');
            end
            y = inline(func);
            plot(t,y(t));
    end

Upvotes: 3

Views: 29009

Answers (2)

macduff
macduff

Reputation: 4685

Looks to me that your problem is what sawtooth returns the mathworks documentation says that:

sawtooth(t,width) generates a modified triangle wave where width, a scalar parameter between 0 and 1, determines the point between 0 and 2π at which the maximum occurs. The function increases from -1 to 1 on the interval 0 to 2πwidth, then decreases linearly from 1 to -1 on the interval 2πwidth to 2π. Thus a parameter of 0.5 specifies a standard triangle wave, symmetric about time instant π with peak-to-peak amplitude of 1. sawtooth(t,1) is equivalent to sawtooth(t).

So I'm guessing that's part of your problem.

After you responded I looked into it some more. Looks to me like it's the quad function; not very accurate! I recast the problem like this:

    function [ a0,am,bm ] = sotest( t, numTerms )
      bm = zeros(1,numTerms);
      am = zeros(1,numTerms);
      % 2L = 1
      L = 0.5;
      for ii = 1:numTerms
        am(ii) = (1/L)*quadl(@(x) aCos(x,ii,L),0,2*L);
        bm(ii) = (1/L)*quadl(@(x) aSin(x,ii,L),0,2*L);
      end
      ii = 0;
      a0 = (1/L)*trapz( t, t.*cos((ii*pi*t)/L) );  
      % now let's test it
      y = ones(size(t))*(a0/2);
      for ii=1:numTerms
        y = y + am(ii)*cos(ii*2*pi*t);
        y = y + bm(ii)*sin(ii*2*pi*t);
      end
      figure; plot( t, y);
    end

    function a = aCos(t,n,L)
      a = t.*cos((n*pi*t)/L);
    end

    function b = aSin(t,n,L)
      b = t.*sin((n*pi*t)/L);
    end

And then I called it like:

[ a0,am,bm ] = sotest( t, 100 );

and I got:

Fourier Series

Sweetness!!!

All I really changed was from quad to quadl. I figured that out by using trapz which worked great until the time vector I was using didn't have enough resolution, which led me to believe it was a numerical issue rather than something fundamental. Hope this helps!

Upvotes: 3

roadrunner66
roadrunner66

Reputation: 7941

To troubleshoot your code I would plot the functions you are using and investigate, how the quad function samples them. You might be undersampling them, so make sure your minimum step size is smaller than the period of the function by at least factor 10.

I would suggest using the FFTs that are built-in to Matlab. Not only is the FFT the most efficient method to compute a spectrum (it is n*log(n) dependent on the length n of the array, whereas the integral in n^2 dependent), it will also give you automatically the frequency points that are supported by your (equally spaced) time data. If you compute the integral yourself (might be needed if datapoints are not equally spaced), you might calculate frequency data that are not resolved (closer spacing than 1/over the spacing in time, i.e. beyond the 'Fourier limit').

Upvotes: 0

Related Questions