Contents

%Tutorial on Fourier Transforms:

    N=1024; %number of points in the time series

    OutFile1='PureToneInteger.dat';
    OutFile2='PureToneNonInteger.dat';
    OutFile3='Aliasing.dat';
    OutFile4='PureTone_Noise.dat';
    OutFile5='Rectangular.dat';
    OutFile6='IrregularSampling.dat';
    OutFile7='Windowed.dat';

    t=0:N-1;    %time bin index array
    f=0:N-1;    %frequency bin index array

Fourier transform of a pure cosine wave (integer frequency)

    figure(1)
    f0=32;          %frequency of the cosine wave
    y=cos(2*pi*t*f0/N);
    z=fft(y);
    Y=z.*conj(z)/(N*N);   %power spectrum

    fid=fopen(OutFile1,'w');
    fprintf(fid,'%f\n',y);
    fclose(fid);

    subplot(3,1,1)
        plot(t,y)
        ylim([-2 2])
        xlim([0 N])
        xlabel('time')
        ylabel('y(t)')
        title(['frequency, f_0= ' num2str(f0)])
        set(gca,'FontSize',16,'LineWidth',2)
    subplot(3,1,2)
        plot(f,Y)
        ylim([0 0.5])
        xlim([0 N])
        ylabel('Power Spectrum, P(f)')
        xlabel('Frequency (Positive & Negative)')
    subplot(3,1,3)
        semilogy(f(1:N/2),Y(1:N/2),'.')
        ylim([1e-12 1])
        xlim([0 N/2])
        legend('note log scale')
        ylabel('log10(P(f))')
        xlabel('Positive Frequency')

Fourier transform of a pure cosine wave (non-integer frequency)

    figure(2)
    f0=100.5;
    y=cos(2*pi*t*f0/N);
    z=fft(y);
    Y=z.*conj(z)/(N*N/4);

    fid=fopen(OutFile2,'w');
    fprintf(fid,'%f\n',y);
    fclose(fid);

    subplot(3,1,1)
        plot(t,y)
        title(['f_0=' num2str(100.5)])
        ylim([-2 2])
        xlim([0 N])
        xlabel('time')
        ylabel('y(t)')
        set(gca,'FontSize',16,'LineWidth',2)
    subplot(3,1,2)
        plot(f(1:N/2),Y(1:N/2))
        xlim([0 N/2])
        ylabel('P(f)')
        xlabel('Positive Frequency')
    subplot(3,1,3)
        plot(f(80:120),Y(80:120),'+','MarkerSize',6)
        xlabel('Positive Frequency')
        ylabel('P(f)')
        xlim([80 120])

Aliasing (what happens if your input time series is high frequency)

    figure(3)
        af0=100;
        f0=N-af0;    %any frequency f_0>N/2 is aliased
        y=cos(2*pi*t*f0/N);
        z=fft(y);
        Y=z.*conj(z)/(N*N/4);

    fid=fopen(OutFile3,'w');
    fprintf(fid,'%f\n',y);
    fclose(fid);

    subplot(2,1,1)
        plot(t,y)
        title(['f_0=' num2str(f0) ' which is' num2str(N) '-' num2str(af0)])
        ylabel('y(t)=cos(2\pi f_0t/N)','FontSize',12)
        xlabel('time')
        xlim([0 N])

    subplot(2,1,2)
        plot(f(1:N/2),Y(1:N/2))
        title(['Frequency is now aliased to freq bin ' num2str(af0)])
        ylabel('Power')
        xlabel('Positive Frequencies')
        xlim([0 N/2])

We now add Gaussian noise to an underlying cosine

    figure(4)
    f0=64;  %frequency of cosine
    a=1;    %rms of Gaussian noise  (rms)

    y=cos(2*pi*t*f0/N);  %amplitude of cosine
    n=a*randn(1,N);  %mean=0, rms=1, Gaussian noise
    yn=y+n;

    fid=fopen(OutFile4,'w');    %write out the signal
    fprintf(fid,'%f\n',yn);
    fclose(fid);

        subplot(3,1,1)
            plot(t,y)
            ylim([-2 2])
            title(['pure tone, f_0=' num2str(f0)])
            xlabel('time')
            ylabel('y(t)')
            xlim([0 N])

        subplot(3,1,2)
            plot(t,yn)
            ylim([-2 2])
            xlim([0 N/2])
            ylabel('y(t)+n(t)')
            title('added unit variance Gaussian noise: ratty signal')
        subplot(3,1,3)
            z=fft(yn);
            Yn=z.*conj(z)/(N*N/4);
            semilogy(f(1:N/2),Yn(1:N/2))
            title('switched to log. Notice the noise floor')
            xlim([0 N/2])
            ylim([1e-3 1])
            ylabel('P(f)')

Periodic signals which are not pure sine (or cosine)

Square wave
    figure(5)
    f0=16.001;
    y=GenerateSignal(N,f0,0);  %my code to generate rect signal

    fid=fopen(OutFile5,'w');    %write out the signal
    fprintf(fid,'%f\n',y);
    fclose(fid);

        subplot(3,1,1)
            plot(t,y);
            ylim([-2 2])
            ylabel('amplitude')
            xlabel('time')
        subplot(3,1,2)
            z=fft(y);
            Y=z.*conj(z)/(N*N/4)+eps;
            plot(f(1:N/2),Y(1:N/2))
            legend('Note the presence of odd harmonics')
            ylabel('P(f)')
            xlabel('f')
        subplot(3,1,3)
            semilogy(f(1:N/2),Y(1:N/2))
            ylim([1e-3 1])
            legend('Note log scale')
            ylabel('P(f)')
            xlabel('f')

Irregular sampling

    figure(6)
        f0=16;
        cutoff=0.8;
        x=cos(2*pi*t*f0/N);
        a=(sign(rand(1,N)-cutoff)+1)/2;  %select only 1/10 data points
        y=a.*x;

        fid=fopen(OutFile6,'w');    %write out the signal
        fprintf(fid,'%f\n',y);
        fclose(fid);

        z=fft(y);
        Y=z.*conj(z)/(N*N/4)+eps;
        subplot(3,1,1)
            plot(t,x,'.')
            xlim([1 N])
            xlabel('time')
            ylabel('amplitude')
            title('pure tone')
        subplot(3,1,2)
            plot(t,y,'.')
            nselect=sum(y~=0);
            legend(['# points=' num2str(nselect)]);
            xlabel('t')
            ylabel('amplitude')
            xlim([0 N])
            title('sparse sampling (1 in 5)')
        subplot(3,1,3)
            semilogy(f(1:N/2),Y(1:N/2))
            xlabel('f')
            ylabel('P(f)')
            xlim([0 N/2])
            ylim([1e-4 1])
            legend('note log scale')

Windowing a sinusoidal signal

    figure(7)
        f0=70;
        F0=4;
        y=WindowedSineWave(N,f0,F0);

        fid=fopen(OutFile7,'w');    %write out the signal
        fprintf(fid,'%f\n',y);
        fclose(fid);

        z=fft(y);
        Y=z.*conj(z)/(N*N/4)+eps;
        subplot(3,1,1)
            plot(t,y)
            xlabel('t')
            ylabel('amplitude')
            title(['f_0=' num2str(f0) ' F_0=' num2str(F0)]);
        subplot(3,1,2)
            plot(f(1:N/2),Y(1:N/2))
            xlabel('f')
            ylabel('P(f)')
        subplot(3,1,3)
            plot(f(1:round(2*f0)),Y(1:round(2*f0)),'+-')
            xlabel('f')
            ylabel('P(f)')
            title('Side lobes from window function')

GenerateSignal

%Here you will find all the routines used in this program

% <include>GenerateSignal.m</include>

% <include>WindowedSineWave</include>