% ps04_1.m is the m-file answer for Problem set 5 question 1
% created by T. Kenna 981015 1148
% modif'd by D. Glover 2000:12:01
% modif'd by D. Glover 2002:10:28
% modif'd by D. Glover 2004:12:20 
% 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('After loading freqanal.dat, we define the measurement');
disp('frequency, fm, and then we create 3 specplots in 1 figure');
disp('using the following commands:');
disp(' ');
disp('fm=1/.005;');
disp('subplot(311)');
disp('specplot(spectrum(freqanal,256),fm)');
disp('subplot(312)');
disp('specplot(spectrum(freqanal,128),fm)');
disp('subplot(313)');
disp('specplot(spectrum(freqanal,64),fm)');
disp(' ');
disp('Let''s take a look at the plot. <return>');
disp(' ');
pause
 
load freqanal.dat
fm=1/.005;
figure(1);clf
subplot(311)
specplot(spectrum(freqanal,256),fm)
text(5,1e4,'Bin = 256')
subplot(312)
specplot(spectrum(freqanal,128),fm)
text(5,1e4,'Bin = 128')
subplot(313)
specplot(spectrum(freqanal,64),fm)
text(5,1e4,'Bin = 64')
pause(life)
% delete(figure(1))

% clc
disp('Notice how as you decrease the number of bins, you get a');
disp('spreading out of the frequency resolution (your ability to');
disp('distinguish two very close together peaks) and a smoothing');
disp('out of the baseline (and hence our ability to "see" a peak');
disp('against the background variability). <return>');
disp(' ');
pause

% clc
disp('Now, to determine the positions of the peaks, you need to');
disp('create a frequency vector matching the frequency bins.');
disp('What you generate is a range of frequencies from 0 to the');
disp('nyquist frequency (fn), which is 1/(2*deltaT) = 100 Hz. Now');
disp('you created (in the first case) 256 bins, which go from');
disp('-fn to +fn, but the only frequencies of interest are');
disp('positive, and "spectrum" gives you Nbins/2 + 1 values');
disp('back, ordered in frequency from 0 to fn in steps of');
disp('fn/128. Thus we create a frequency vector to match the');
disp('bins, sort them in descending order of power (that''s');
disp('why the negative sign), save the indices (ind) into');
disp('the sorted vector and use this index vector to list');
disp('the top ten bins. Note: Column 1 is Frequency, Column 2');
disp('is power spectral density, and column 3 is uncertanity'); 
disp(' ');
disp('fn=fm/2;');
disp('p=spectrum(freqanal,256);');
disp('f=[0:fn/128:fn]'';');
disp('[fs,ind]=sort(-p(:,1));');
disp('Top_Ten_Bins = [f(ind(1:10)) p(ind(1:10),:)]');
% disp('Top_Ten_Bins = flipud(sort(Top_Ten_Bins))');
disp(' ');
disp('<return>');
pause

% clc
fn=fm/2;
f=[0:fn/128:fn]';
p=spectrum(freqanal,256);
[fs,ind]=sort(-p(:,1));
Top_Ten_Bins=[f(ind(1:10)) p(ind(1:10),:)];
disp('     Frequency       Power   Uncertainty (PSD)');
% Top_Ten_Bins = flipud(sort(Top_Ten_Bins))
Top_Ten_Bins
disp('Now the three peak are centered on 85 Hz, 54 Hz, and 16 Hz.');
disp('An alternate approach to finding the peak is to simply list');
disp('the values with a command like "[f p]" and then scroll');
disp('through the list to find the peaks.'); 
disp(' ');
pause

% clc
disp('Next, we add the frequencies 110 hz, 160 hz and 146 hz and plot');
disp('up the spectra. The results are shown below. <return>');
disp(' ');
pause
 
figure(2);clf
t=0:.005:10;
y=freqanal';
f=110;
y=y + 3*sin(2*pi*f*t);
specplot(spectrum(y,256),fm)
text(60,7e2,'The aliased 110 hz peak')
text(60,4e2,'shows up at 90 hz')
text(5,5e-3,'f=110')
x=1e-3:100:1e3;
y=89.5*ones(size(x));
h=line(y,x);
set(h,'color','red','Linewidth',2)  
pause(life)
% delete(figure(2))

% clc
disp('Notice that the frequency is above the nyquist');
disp('frequency, so it is reflected back toward zero,');
disp('and appears at 90 Hz (since it was 10 Hz above');
disp('the nyquist frequency, it shows up at 10 Hz');
disp('below the nyquist frequency). This can be');
disp('represented by fa = 2fn - f where "fa" is the');
disp('apparent frequency, "f" is the actual');
disp('frequency, and "fn" is the nyquist frequency.');
disp('Next, let''s look at the 160 frequency. <return>');
disp(' ');
pause 

figure(3);clf
t=0:.005:10;
y=freqanal';
f=160;
y=y + 3*sin(2*pi*f*t);
specplot(spectrum(y,256),fm)
text(42,7e2,'The aliased 160 hz peak')
text(42,4e2,'shows up at 40 hz')
text(5,5e-3,'f=160')
x=1e-3:100:1e3;
y=40*ones(size(x));
h=line(y,x);
set(h,'color','red','Linewidth',2)  
pause(life)
% delete(figure(3))

% clc
disp('Now the 160 Hz signal is even further above the');
disp('nyquist frequency, so it is reflected back by its');
disp('distance above the nyquist limit and shows up at');
disp('40 Hz now: (since it is 60 Hz above the nyquist');
disp('frequency, it will now show up at 2fn -f = 40 Hz).');
disp('Finally. let''s look at the 146 frequency. <return>');
disp(' ');
pause 

figure(4);clf
t=0:.005:10;
y=freqanal';
f=146;
y=y + 3*sin(2*pi*f*t); 
specplot(spectrum(y,256),fm)
text(22,7e2,'The aliased 146 hz peak')
text(22,4e2,'shows up at 54 hz and')
text(22,2e2,'combines destructively')
text(22,110,'with the existing 54 hz')
text(22,70,'peak')
text(5,5e-3,'f=146')
x=1e-3:100:1e3;
y=54*ones(size(x));
h=line(y,x);
set(h,'color','red','Linewidth',2)  
pause(life)
% delete(figure(4))

% clc
disp('Now the situation becomes a little tricky. The 146');
disp('Hz frequency is reflected back onto the existing 54');
disp('Hz peak ( 2fn - 146 = 54). Hmm, something''s');
disp('interesting here: the 54 Hz peak actually decreased');
disp('in size! You''d think that when you add another peak');
disp('to it that it would increase, not decrease. Well');
disp('actually, what is happening is that the two frequencies');
disp('are interfering, and that the "amplitude of their sums');
disp('is not equal to the sum of their amplitudes". You can');
disp('show this by adding two sine waves of different');
disp('frequencies. <return>'); 
disp(' ');
pause
 
% clc
disp('End of Problem 1');