% eda06_04
% amplitude spectral density of the Neuse River hydrograph
% (not normalized - see eda06_13 for properly normalized version)

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;
M=N;
d=d(1:N);

% set up time
tmin=0.0;
tmax = tmin+Dt*(N-1);
t = tmin + Dt*[0:N-1]';

% nyquist frequencies
Nf=N/2+1;
fmax = 1/(2*Dt);
Df = fmax/(N/2);
f = Df*[0:Nf-1]';
Nw=Nf;
wmax = 2*pi*fmax;
Dw = wmax/(N/2);
w = Dw*[0:Nw-1]';

% set up G
G=zeros(N,M);

% zero frequency column
G(:,1)=1;

% interior M/2-1 columns
for i = [1:M/2-1]
    j = 2*i;
    k = j+1;
    G(:,j)=cos(w(i+1).*t);
    G(:,k)=sin(w(i+1).*t);
end

% nyquist column
G(:,M)=cos(w(Nw).*t);
    
% solve least squares problem
gtgi = 2*ones(M,1)/N;
gtgi(1)=1/N;
gtgi(M)=1/N;
mest = gtgi .* (G'*d);

% check results
dpre=G*mest;
E=(d-dpre)'*(d-dpre);
E

% compute spectral density

% zero frequency
s=zeros(Nw,1);
s(:,1)=sqrt(mest(1)^2);

% interior points
for i = [1:M/2-1]
    j = 2*i;
    k = j+1;
    s(i+1) = sqrt(mest(j)^2 + mest(k)^2);
end

% Nyquist frequency
s(Nw) = sqrt(mest(M)^2);

% plot spectral density
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 ] );
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 ] );
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]);
hold on;
plot(1./f(2:Nf),s(2:Nf),'k-','LineWidth',2);
xlabel('period, days');
ylabel('spectrum');



