# Computing the power spectrum via matlab

1. Sep 19, 2009

### oscarmike

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
t = timescale (:,1);
p = timescale (:,2);

clf, axis % clear graph and reset plot axis
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

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?