Contents

read pulsar data

    infile='psr_data/medium_psr.dat';
    fid=fopen(infile);

    yraw=fread(fid,'single');
    fclose(fid);
    n=numel(yraw);
    t=1:n;

"Condition" the data

%astronomical receivers have small gain changes
%there are also changes in the background contribution
%NSMOOTH is the smoothing timescale. 10000 corresponds
%to 10 seconds (in this case).

    NSMOOTH=10000;
    ysmooth=smooth(yraw,NSMOOTH);
    y=yraw-ysmooth; %subtract the smoothed data
    yvar=var(y);    %variance of the input time series
    figure(1)
    plot(t,yraw,'.')
    hold on
    plot(t,ysmooth,'-k')
    hold off

    xlabel('Bin')
    ylabel('y')
    title('Raw Data & Smoothed Curve')
    set(gca,'FontSize',16,'LineWidth',2)

Fourier Transform

FFT is optimized when number of data points is a product of prime numbers, especially 2^n

% We increase the size of the array to accommodate the
% next power of 2
    N=2^(nextpow2(n));

    Y=fft(y,N);
    P=Y.*conj(Y)/(yvar*N);   %Normalization
    P(N/2:N)=[];            %No need to retain -ive freqs
    P(1:100)=P(101);        %lots of junk at low freq
                            %suppress it

Plot the data

    figure(2)
    semilogy(P,'.')
    hold on
    xl=4000;
    xlim([1 xl])
    ylim([0.1 3E2])

    xlabel('Frequency bin')
    ylabel('Power Spectrum')
    set(gca,'FontSize',16,'LineWidth',2)

Mark harmonics

    f0=287.09;  %I determined the fundamental by looking
                %at the power spectrum
    mharm=11;
    ind=round(f0*(1:mharm)+1);
    plot(ind,P(ind),'ok','MarkerSize',10)
    title('Harmonics circled')
    hold off

Measure the bin location of the peaks

    figure(3)
    mharm=10;
   [pharm,binharm]=FindHarmonics(y,round(f0),mharm);
   a=binharm./(1:mharm);
   plot(1:mharm,a,'+','MarkerSize',12)
   xlabel('harmonic number')
   ylabel('frequency')
   xlim([0 11])