% ps04_2.m is the m-file answer for Problem set 5 question 2
% created by T. Kenna 981015 1150
% Modif'd by D. Glover 2000:12:01
% Modif'd by D. Glover 2000:12:08 to add to avg Nyquist frequency
% Modif'd by D. Glover 2002:10:28
% Modif'd by D. Glover 2008:10:20 mv to PS04 from PS05
%

clear all
close all

global G
global g
global n

life=0;

% clc
disp('First, we load the data and  calculate the average nyquist');
disp('frequency using the following commands (according to Press et al.)');
disp(' ');
disp('load tu.dat');
disp('load fu.dat');
disp('dif=max(tu)-min(tu)');
disp('Average_delta_t=dif/size(tu,1)');
disp('Average_Nyquist_Freq=inv(2*Average_delta_t)'); 
disp(' ');
disp('<return>');
pause

load tu.dat
load fu.dat
dif=max(tu)-min(tu)
Average_delta_t=dif/size(tu,1)
Average_Nyquist_Freq=inv(2*Average_delta_t) 

disp('Or we could take a slightly different approach by calculating');
disp('the average delta-t between samples. This has the advantage of');
disp('weighting your sampling rate if you have a bunch of data');
disp('collected at short delta-t''s PLUS a few more widely spaced');
disp('delta-t''s. In this case we would calculate the average Nyquist');
disp('frequency in the following manner:');
disp(' ');
disp('diff_t=diff(tu)');
disp('Average_delta_t=mean(diff_t)');
disp('Average_Nyquist_Freq=inv(2*Average_delta_t)'); 
disp(' ');
disp('<return>');
pause

diff_t=diff(tu);
Average_delta_t=mean(diff_t)
Average_Nyquist_Freq=inv(2*Average_delta_t)

disp('As you can see, in this case, they are very close.');
disp(' ');
disp('<return>');
pause

disp('If we plot the original time vs observation, we can');
disp('see the irregularity of the interval.');
disp(' ');
pause

figure(1);clf
plot(tu,fu,'r*')                                                  
title('Input Data')
xlabel('Time (sec)')
ylabel('Observation')
pause(life)
% delete(figure(1))

% clc
disp('If we make a histogram plot of the log10(separation), we observe');
disp('that the bulk of the observations are separated by an interval');
disp('of around .01 however there are almost an equal number of');
disp('observations that are close to .003 as well as a smattering of');
disp('points with even finer spacing. <return>');
disp(' ');
pause

figure(2);clf
[N,X]=hist(log10(diff(tu)),20);
bar(X,N,'b');                                                    
title('Intervals')                    
xlabel('Log10(deltat)')
pause(life)
% delete(figure(2))

% clc
disp('We now use the Lomb method to create a normalized PSD, and');
disp('probability whether the peaks are random or not; we plot both');
disp('against a frequency vector that goes from zero to approximately');
disp('3 times the average nyquist frequency (this decision is somewhat');
disp('arbitrary, you can go out as far as you want, you may observe peaks,');
disp('but their probabilities will be close to one. (if Prob=1; there is a');
disp('100% probability that the peak is random). <return>');
disp(' ');
pause
 
figure(3);clf
subplot(211)
freq=1:1:150;   
[Pn, Prob]=lomb(tu,fu,freq);
plot(freq,Pn,'r')                                            
title('Lomb Periodogram') 
xlabel('Frequency (Hz)')
ylabel('Normalized PSD')
subplot(212)
plot(freq,Prob,'b')       
ylabel('Probability')
xlabel('Frequency (Hz)')
pause(life)
% delete(figure(3))

% clc
disp('We can create a matrix of Frequency, PSD and probablility and');
disp('then sort them according to some criteria; in this case I only');
disp('want the frequencies that have 80% or less probability of being');
disp('random. <return>');
pause

% clc       
Table=[freq' Pn' Prob'];        
Information_Table=sortrows(Table(Table(:,3)<=.8,1:3))
disp('The dominant peaks returned are 12, 42, 47, and 70.  The probability');
disp('that each is random is 3.09%, 0%, 73.81%, and  1.16% respectively.');
disp('Because there is only a 26.19% chance that the 47 Hz peak is real, it');
disp('should not be considered. If we compare our results to the average');
disp('nyquist frequency of 50.4154, we have one frequency(12) which is');
disp('below the average nyquist value, another frequency(42) which is near');
disp('the average nyquist value and the third frequency (70) which exceeds');
disp('the average nyquist value.  We see the benefit of an irregularly');
disp('spaced data set in that, had we spaced the measurements at regular');
disp('points, we may have been able to calculate the first two frequencies');
disp('however we would have been plagued by aliasing with the 70 Hz');
disp('frequency. <return>');
disp(' ');
pause

% clc
disp('End of Problem 2');