1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Computing the power spectrum via matlab

  1. Sep 19, 2009 #1
    1. The problem statement, all variables and given/known data
    I need to firstly fourier transform some data and then compute and plot the power spectrum. The topic is modelling magnetic anomalies. My code is below:

    Code (Text):
    % Matlab script to load and plot bathymetry and magnetics
    % profile from southeast Indian Ridge
    % The first column of the file is age (in my)
    % The second column of the file is normal (+1) or
    % reversed (-1) field polarity
    load timescale.dat
    t = timescale (:,1);
    p = timescale (:,2);

    clf, axis % clear graph and reset plot axis
    % load timescale.dat
    t =timescale(:,1);
    p =timescale(:,2);
    %
    % plot timescale
    %
    subplot(2,1,1)
    plot(t,p)
    xlabel(' Million Years Before Present')
    ylabel(' Magnetic Polarity')
    title (' Magnetic Timescale for 10 Myr')
    axis([0.,10.,-1.5,1.5]);
    pause

    load sei_mag.dat
    load sei_topo.dat
    lat =sei_mag(:,2);
    lat_topo =sei_topo(:,2);
    magobs =sei_mag(:,3);
    depth =sei_topo(:,3);
    subplot(2,1,1)
    plot(lat,magobs)
    ylabel('Magnetic Anomaly (n tesla)')
    xlabel('Latitude (deg)') %Not sure if this is correct.
    title (' Magnetic Data')

    subplot(2,1,2)
    plot(lat_topo, depth)
    xlabel('Latitude (deg)')
    ylabel('Depth (m)')
    title (' Topography Data')
    pause


    [ndat,mdat]=size(t);
    mt=-fliplr(t(2:ndat)')';
    mp= fliplr(p(2:ndat)')';
    time=[mt;t];
    polarity=[mp;p];

    %Question 6
    [ndat,mdat] = size(time);
    nmid = round(ndat/2);
    %
    % time is in increments of 10k,
    % so +/- 10My will be a 2048 long vector
    %
    ntime=2048;
    ntime2=ntime/2;
    stime = time((nmid-ntime2+1):(nmid+ntime2),1);
    spolr = polarity((nmid-ntime2+1):(nmid+ntime2),1);
    plot(stime',spolr')
    axis([-10,10,-1.2,1.2]);
    xlabel(' Million Years Before Present')
    ylabel(' Magnetic Polarity')
    pause
    axis;


    % halfrate v = 30000 m/my; dt = 20 my
    % L =v*dt is in meters, a = 3000 m
    %
    v = 30000;
    dt = 20;
    L = v*dt;
    nx = 2048;
    nx2 = nx/2;
    dx = L/(nx-1);
    k = -nx2:(nx2-1);
    k = k' ./ L;

    %Question 8 - Power Spectrum
    % Fourier transform data data into the frequency domain


    spolrd=detrend(spolr);
    stimed=detrend(stime);

    spolr1=fft(spolrd,nx);
    stime1=fft(stimed, nx);

    YTP=(2*dt/nx)*spolr1 .* conj(spolr1);   %raw power spectrum
    semilogx(k,abs(YTP(1:nx2)))
    grid on
    title('Cross-spectrum of magnetic data')
    xlabel('TBA')
    ylabel('Cross-spectral power')

     

    But the code that I am having trouble with is:

    Code (Text):
    %Question 8 - Power Spectrum
    % Fourier transform data data into the frequency domain


    spolrd=detrend(spolr);
    stimed=detrend(stime);

    spolr1=fft(spolrd,nx);
    stime1=fft(stimed, nx);

    YTP=(2*dt/nx)*spolr1 .* conj(spolr1);   %raw power spectrum
    semilogx(k,abs(YTP(1:nx2)))
    grid on
    title('Cross-spectrum of magnetic data')
    xlabel('TBA')
    ylabel('Cross-spectral power')
     
    I can't get it to plot as it seems the vectors are different lengths, does that mean I have done the fourier transform wrong?
     
  2. jcsd
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?