% eda06_05
% amplitude spectral density of the Neuse River hydrograph
% using fft command
% (unnormalized spectral density - see eda06_13 for proper normalization)

clear all;

D = load('neuse.txt');
Dt = D(2,1)-D(1,1);
d = D(:,2);
N=length(d);

% round off to even number of points
N=floor(N/2)*2;
d=d(1:N);

% generic time/frequency set up
M=N;
tmax=Dt*(N-1); 
t=Dt*[0:N-1]'; 
fmax=1/(2.0*Dt);
df=fmax/(N/2);
f=df*[0:N/2,-N/2+1:-1]'; 
Nf=N/2+1;
dw=2*pi*df;
w=dw*[0:N/2,-N/2+1:-1]'; 
Nw=Nf;

% compute fourier coefficients
mest = fft(d);

% compute amplitude spectral density
s=abs(mest(1:Nw));

% plot spectrum
figure(1)
clf;

% plot frequency
subplot(4,1,1);
set(gca,'LineWidth',2);
hold on;
plot(f(2:Nf),s(2:Nf),'k-','LineWidth',2);
axis( [0, 0.5, 0, 2000*N/2 ] );
xlabel('frequency, cycles per day');
ylabel('spectrum');

% plot log in frequency
subplot(4,1,2);
lw = 2.2;
set(gca,'LineWidth',lw,'Xscale','log','Yscale','lin');  hold on;
plot(f(2:Nf),s(2:Nf),'k-','LineWidth',2);
axis( [.001, 0.5, 0,2000*N/2 ] );
xlabel('frequency, cycles per day');
ylabel('spectrum');

% plot period
subplot(4,1,3);
set(gca,'LineWidth',2);
hold on;
axis( [0, 1000, 0, max(s(2:Nf)) ] );
plot(1./f(2:Nf),s(2:Nf),'k-','LineWidth',2);
xlabel('period, days');
ylabel('spectrum');

% log plot in period
subplot(4,1,4);
lw = 2.2;
set(gca,'LineWidth',lw,'Xscale','log','Yscale','lin');
hold on;
axis([10,1000,0,2000*N/2]);
hold on;
plot(1./f(2:Nf),s(2:Nf),'k-','LineWidth',2);
xlabel('period, days');
ylabel('spectrum');

% compute original data from fourier coefficients
dnew = ifft(mest);
E = (d-dnew)'*(d-dnew);
fprintf('Error in d->mest->dnew calculation: %f\n', E );