% Program to Calculate the Cross Power Spectrum of Two Channels of Data

% Change Directories
cd 'C:\Data'

% Load Recorded Data File (Model.txt)
load Model.txt 

% Define First Column of File as Time Array
time = Model(:,1); 

% Define Second Column of File as Channel 1 Acceleration Data Array
ch1 = Model(:,2); 

% Define Third Column of File as Channel 2 Acceleration Data Array
ch2 = Model(:,3); 

Remove Linear Offset from Time History by Detrending
ch1 = detrend(ch1);
ch2 = detrend(ch2);

% Compute the Fast Fourier Transform (FFT) of Channel 1 Data
ch1_fft = fft(ch1); 

% Determine the Length of the FFT of Channel 1
N = length(ch1_fft); 

% Calculate delta t by Subtracting Point 2 from Point 1 in the Time Array
deltat=time(2)-time(1) 

% Determine the Number of Time Steps (nts)
nts=length(ch1) 

% Calculate the Frequency Increment of the FFT
deltaf=1/nts/deltat 

% Calculate the Nyquist Frequnecy
nyq=1/deltat/2 

% Set up the Frequency Array Which Ranges From 0 to the Nyquist Frequency in Delta f Increments
f=0:deltaf:nyq-deltaf; 

% Calculate the FFT of Channel 2 Data
ch2_fft = fft(ch2); 

% Calculate the Complex Conjugate of the FFT of Channel 2
ch2_fft_conj = conj(ch2_fft); 

% Calculate the Number of Time Steps Sqaured
nts2=nts^2; 

% Start Calculating the Cross Power Spectrum
cross=ch1_fft.*ch2_fft_conj; 
CrossSpectrum=cross./nts2; %Finish Calculating the Cross Power Spectrum


% Convert from a Double-Sided Spectrum to a Single-Sided Spectrum
CrossSpectrum=CrossSpectrum(1:N/2); %Discard all points from i = N/2 + 1 through N
CrossSpectrum=CrossSpectrum.*2; %Multiply all points from i = 1 to N/2 by 2
CrossSpectrum(1)=CrossSpectrum(1)./2; %The DC point should not be multiplied by 2

% Compute the Cross Power Spectrum Magnitude
mag_CrossSpectrum = abs(CrossSpectrum(1:N/2)); 

% Compute the Cross Power Spectrum Phase Angle
theta_CrossSpectrum = angle(CrossSpectrum(1:N/2)); 

% Plot the Cross Power Spectrum Magnitude versus frequency
plot(f,mag_CrossSpectrum)
xlabel('Frequency (Hz)') % X-axis Label
ylabel('g (RMS)') % Y-axis Label
title('Cross Power Spectrum Magnitude Versus Frequency') % Plot Title

% Plot the Cross Power Spectrum Phase Angle versus frequency
plot(f,theta_CrossSpectrum)
xlabel('Frequency (Hz)') % X-axis Label
ylabel('Phase Angle (Radians)') % Y-axis Label
title('Cross Power Spectrum Phase Angle Versus Frequency') % Plot Title