% frequencyanalysis_short.m   WCR 20080418
% Frequency analysis in Matlab - short version.
% Reads data in from file.  Plots the data and the detrended, windowed data.
% Computes the Fourier transform, plots the scaled amplitudes and phases.
% Computes and plots the one-sided power spectrum.
% Amplitudes are "two-sided" i.e. 0.707xLabview's "one-sided" amplitudes. 
% One-sided power spectrum matches Labview's one-sided power spectrum. 

filespec=input('Enter file to read in, should have 1 column of data: ','s');
y1=load(filespec);          % y1 = data from file
N=length(y1);               % how many data points
fs=input('Enter sampling rate in Hz: ');    % get sampling rate from user
dt=1/fs;                    % sampling interval (in seconds)
t=dt*[0:N-1];               % times to use as x-values for plotting
y=detrend(y1).*hamming(N);  % y = detrended, Hamming-windowed data
plot(t,y1,'b',t,y,'r');     % plot raw & detrended data (blue & red)
grid on; legend('raw data','detrended windowed data');   % plot grid & legend
xlabel('Time (s)'); ylabel('Y axis (units)');   % add axis labels
yfft=fft(y)/N;              % yfft = scaled Fourier transform of detrended windowed data
f=[0:fs/N:fs/2];            % frequencies to use as x values for plotting
yfftmag=abs(yfft(1:(N/2+1)));   % amplitudes of scaled Fourier transform (FT)
yfftphs=angle(yfft(1:(N/2+1))); % phase angles (in radians) of the FT
figure; subplot(2,1,1);     % prepare figure with 2 panels for plotting
plot(f,yfftmag);            % plot amplitude vs. frequency
grid on; xlabel('Frequency (Hz)'); ylabel('Amplitude (units)'); %prettify
subplot(2,1,2);             % prepare second panel of figure for plotting
plot(f,yfftphs);            % plot phase  vs. frequency
grid on; xlabel('Frequency (Hz)'); ylabel('Phase (rad)'); % prettify
                            % next: compute 1-sided power spectrum
yps1=[(yfft(1))^2; 2*abs(yfft(2:N/2)).^2; abs(yfft(N/2+1)).^2];
figure;                     % prepare figure for next plot
plot(f,yps1);               % plot 1-sided power spectrum vs frequency
grid on; xlabel('Frequency (Hz)'); ylabel('Power (units^2)'); % prettify
