Matlab+Code

**Question 2.1**
code format="matlab" %% Specify 'flag' to plot figures! clear;clc;close; flag = 1; % 1: plots ch_2.1.1; 2: plots ch_2.1.2; 3: plots ch_2.1.3, etc.

%% 1. grab data: sst, slp sst = nc_varget('data.nc','sst'); slp = nc_varget('data.nc','slp'); ts1 = sst(:,80); ts2 = slp(:,80);

%% 2. FFT T = 240;                  % number of records p = abs(fft(ts1)/T);      % absolute value of the fft, careful with scaling % p(1) = mean(ts1) = 27.8039 --> mean value, zero frequency %__________________________________________________________________________ % 2.1 PSD vs f if flag==1 % 2.1.1 psd with f centers at 0 p = p.^2;              % sum(p(2:end)) == var(ts1,1) =1.1118 --> power = variance f = [-119:120]/T; semilogy(f,cat(1,p(122:end),p(1:121))/(1/T),'-*k',...      'MarkerSize',4,'MarkerEdgeColor','r','LineWidth',0.8); % semilogy(f,fftshift(p)/(1/T),'-*k',...  %    'MarkerSize',4,'MarkerEdgeColor','r','LineWidth',0.8);   % the same, but f=[-120:119]/T; xlabel('frequency (cycle/month)','FontSize',12); ylabel('PSD (^oC^2·month)','FontSize',12); title('Power Spectrum Density','FontSize',12); grid on; set(gca,'FontSize',12); print('-djpeg','ch_2.1.1') elseif flag==2 % 2.1.2 half of the spectrum p = p(2:T/2+1).^2; p = p.*2;              % periods of data been split into 'positive' and 'negative' f = [1:T/2]/T; plot(f,p/(1/T),'-*k','MarkerSize',4,'MarkerEdgeColor','r','LineWidth',0.8); xlabel('frequency (cycle/month)','FontSize',12); ylabel('PSD (^oC^2·month)','FontSize',12); title('Power Spectrum Density','FontSize',12); grid on; set(gca,'FontSize',12); print('-djpeg','ch_2.1.2') elseif flag==3 % 2.1.3 rebin p = p(2:T/2+1).^2; p = p.*2;              % periods of data been split into 'positive' and 'negative' f = [1:T/2]/T; j = 0; for i=1:2:120 j = j+1; pn(j) = (p(i)+p(i+1))/2; fn(j) = (f(i)+f(i+1))/2; end plot(fn,pn/(1/T),'-*k','MarkerSize',4,'MarkerEdgeColor','r','LineWidth',0.8); xlabel('frequency (cycle/month)','FontSize',12); ylabel('PSD (^oC^2·month)','FontSize',12); title('Power Spectrum Density','FontSize',12); grid on; set(gca,'FontSize',12); xlim([fn(1) fn(end)]); print('-djpeg','ch_2.1.3') elseif flag==4 % 2.1.4 smooth3: average over every 3 values p = p(2:T/2+1).^2; p = p.*2;              % periods of data been split into 'positive' and 'negative' for i=2:119 p(i) = (p(i-1)+p(i)+p(i+1))/3;  % smooth3 end p(120) = []; p(1) = []; f = [2:T/2-1]/T; plot(f,p/(1/T),'-*k','MarkerSize',4,'MarkerEdgeColor','r','LineWidth',0.8); xlabel('frequency (cycle/month)','FontSize',12); ylabel('PSD (^oC^2·month)','FontSize',12); title('Power Spectrum Density','FontSize',12); grid on; set(gca,'FontSize',12); print('-djpeg','ch_2.1.4') end

code

**Question 2.2-2.7**
code format="matlab" %% Specify 'flag' to plot figures! clear;clc;close; flag = 2; % 2: plots ch_2.2; 3: plots ch_2.3; 4: plots ch_2.4, etc.

%% 1. grab data: sst, slp sst = nc_varget('data.nc','sst'); slp = nc_varget('data.nc','slp'); ts1 = sst(:,80); ts2 = slp(:,80);

%% 2. FFT T = 240;               % number of records p = abs(fft(ts1)/T);   % absolute value of the fft, careful with scaling % p(1) = mean(ts1) = 27.8039 --> mean value, zero frequency. p = p(2:T/2+1).^2;     % var(ts1,1)=sum(p)-p(end)/2=1.1122-0.0004=1.1118. p = p.*2;              % periods of data been split into 'positive' and 'negative' f = [1:T/2]/T;         % frequency runs from 1/240 to 120/240

if flag==2 %__________________________________________________________________________  % 2.2 acumulative power vs period pc = cumsum(p); pc(end) = pc(end)-p(end)/2;  % pc(end) = 1.1118; equal to the variance var(ts1,1)

plot(1./f,pc/(1/T),'-*k','MarkerSize',4,'MarkerEdgeColor','b','LineWidth',0.8); xlabel('period (month)','FontSize',12); ylabel('PSD (^oC^2·month)','FontSize',12); title('Acumulative Power Spectrum Density','FontSize',12); grid on; set(gca,'FontSize',12); xlim([1/f(end) 1/f(1)]); print('-djpeg','ch_2.2') elseif flag==3 %__________________________________________________________________________  % 2.3 f*power vs log(f) semilogx(f,f'.*p,'-*k','MarkerSize',4,'MarkerEdgeColor','m','LineWidth',0.8); xlabel('f (cycle/month)','FontSize',12); ylabel('power*f (^oC^2·month^-^1)','FontSize',12); title('Power Spectrum','FontSize',12); grid on; set(gca,'FontSize',12); print('-djpeg','ch_2.3') elseif flag==4 %__________________________________________________________________________  % 2.4 f*power vs log(T) semilogx(1./f,f'.*p,'-*k','MarkerSize',4,'MarkerEdgeColor','m','LineWidth',0.8); xlabel('period (month)','FontSize',12); ylabel('power*f (^oC^2·month^-^1)','FontSize',12); title('Power Spectrum','FontSize',12); grid on; set(gca,'FontSize',12); print('-djpeg','ch_2.4') elseif flag==5 %__________________________________________________________________________  % 2.5 log(power) vs log(f) loglog(f,p,'-*k','MarkerSize',4,'MarkerEdgeColor','m','LineWidth',0.8); xlabel('f (cycle/month)','FontSize',12); ylabel('power (^oC^2)','FontSize',12); title('Power Spectrum','FontSize',12); grid on; set(gca,'FontSize',12); print('-djpeg','ch_2.5') elseif flag==6 %__________________________________________________________________________  % 2.6 rebin again, repeat 2.1.3, but average over 3 points instead of 2 j = 0; for i=1:3:117 j = j+1; pn(j) = (p(i)+p(i+1)+p(i+2))/3; fn(j) = (f(i)+f(i+1)+f(i+2))/3; end loglog(fn,pn/(1/T),'-*k','MarkerSize',4,'MarkerEdgeColor','r','LineWidth',0.8); xlabel('frequency (cycle/month)','FontSize',12); ylabel('PSD (^oC^2·month)','FontSize',12); title('Power Spectrum Density','FontSize',12); grid on; set(gca,'FontSize',12); print('-djpeg','ch_2.6') elseif flag==7 %__________________________________________________________________________  % 2.7 overplotting PSD loglog(f,p,'-*k','MarkerSize',4,'MarkerEdgeColor','m','LineWidth',0.8); hold on;

TL = 240*3; tspad = cat(1,ts1*0,ts1-mean(ts1),ts1*0); pL = abs(fft(tspad)/TL); pL = pL(2:TL/2+1).^2; pL = pL.*2; fL = [1:TL/2]/TL; loglog(fL,pL,'-ob','MarkerSize',4,'MarkerFaceColor','g','LineWidth',0.1);

xlabel('f (cycle/month)','FontSize',12); ylabel('power (^oC^2)','FontSize',12); title('Power Spectrum','FontSize',12); grid on; set(gca,'FontSize',12); legend('ts1','tspad'); print('-djpeg','ch_2.7') elseif flag==8 %__________________________________________________________________________  % 2.8 overplotting PSD, revised of 2.7 loglog(f,p,'-*k','MarkerSize',4,'MarkerEdgeColor','m','LineWidth',0.8); hold on;

TL = 240*3; tspad = cat(1,ts1*0,ts1-mean(ts1),ts1*0); pL = abs(fft(tspad)/TL); pL = pL(2:TL/2+1).^2; pL = pL.*2; fL = [1:TL/2]/TL; j = 0; for i=1:3:358 j = j+1; pLn(j) = (pL(i)+pL(i+1)+pL(i+2))/3; fLn(j) = (fL(i)+fL(i+1)+fL(i+2))/3; end loglog(fLn,pLn,'-ob','MarkerSize',4,'MarkerFaceColor','g','LineWidth',0.1);

xlabel('f (cycle/month)','FontSize',12); ylabel('power (^oC^2)','FontSize',12); title('Power Spectrum','FontSize',12); grid on; set(gca,'FontSize',12); legend('ts1','tspad'); print('-djpeg','ch_2.8') %% slope is about 2?! end code

**Question 3**
code format="matlab" clear;clc;close % grab data: sst, slp sst = nc_varget('data.nc','sst'); slp = nc_varget('data.nc','slp'); ts1 = sst(:,80); ts2 = slp(:,80);

N =240; f = [1:N/2]/N;

% power spectrum of red noise [r,lags] = xcorr(ts1,'coeff'); r1 = r(241); T = -1/log(r1); for i=1:size(f,2) P_red(i) = 2*T/(1 + f(i)^2*T^2); end P_red = P_red*(var(ts1)/sum(P_red)); loglog(f,P_red,'-*k','MarkerSize',4,'MarkerEdgeColor','r','LineWidth',0.8) hold on;

% power spectrum of sst p = abs(fft(ts1)/N); p = p(2:N/2+1).^2; p = p.*2; plot(f,p,'-ob','MarkerSize',4,'MarkerfaceColor','g','LineWidth',0.8); plot(f,P_red*4.61,'-m','LineWidth',1.0); % the value 4.16 is from F-test table

% refine figure xlabel('frequency','FontSize',12); ylabel('power','FontSize',12); title('Power Spectrum','FontSize',12); grid on; set(gca,'FontSize',12); legend('red noise','sst','99 percent quantiles of F-test'); print('-djpeg','ch_3') code

**Question 4**
code format="matlab" % Specify 'flag' to plot figures! clear;clc;close; flag = 4; % 1: plots ch_4.1; 4: plots ch_4.4; 5: plots ch_4.5, etc.

% grab data: sst, slp sst = nc_varget('data.nc','sst'); slp = nc_varget('data.nc','slp'); ts1 = sst(:,80); ts2 = slp(:,80);

T = 240;

if flag==1 %__________________________________________________________________________  % 4.1 Power spectrum of sst and slp f = [1:T/2]/T;

xhat = abs(fft(ts1)/T); xhat = xhat(2:T/2+1).^2; xhat = xhat.*2;

yhat = abs(fft(ts2)/T); yhat = yhat(2:T/2+1).^2; yhat = yhat.*2;

subplot(1,2,1); loglog(f,xhat,'-+k','MarkerSize',4,'MarkerEdgeColor','r','LineWidth',0.8); hold on; xlabel('frequency (cycle/month)','FontSize',12); ylabel('power (^oC^2)','FontSize',12); title('Power Spectrum of SST','FontSize',12); grid on; axis square; set(gca,'FontSize',12);

subplot(1,2,2); loglog(f,yhat,'-bo','MarkerSize',4,'MarkerFaceColor','g','LineWidth',0.8); xlabel('frequency (cycle/month)','FontSize',12); ylabel('Power (^oC^2)','FontSize',12); title('Power Spectrum of SLP','FontSize',12); grid on; axis square; set(gca,'FontSize',12); print('-djpeg','ch_4.1') elseif flag==4 %__________________________________________________________________________  % 4.2 & 4.3 & 4.4 cross-spectrum xhat = fft(ts1)/T; yhat = fft(ts2)/T; cross = xhat .* conj(yhat); R = real(cross);   % in-phase part P (or R)   I = imag(cross);    % quadrature part Q (or I)

R = cat(1,(R(2:120)+flipud(R(122:end))),R(121)); Rc = cumsum(R);    % Rc(end) = -0.1392 = mean((ts1-mean(ts1)).*(ts2-mean(ts2)))

f = [1:T/2]/T; plot(1./f,Rc/(1/T),'-*k','MarkerSize',4,'MarkerEdgeColor','b','LineWidth',0.8); xlabel('period (month)','FontSize',12); ylabel('PSD (^oC^2·month)','FontSize',12); title('Acumulative Power','FontSize',12); grid on; set(gca,'FontSize',12); xlim([1/f(end) 1/f(1)]); print('-djpeg','ch_4.4') elseif flag==5 %__________________________________________________________________________  % 4.5 squared coherency spectrum f = [1:T/2]/T;

xhat = fft(ts1)/T; yhat = fft(ts2)/T; cross = xhat .* conj(yhat); R = real(cross);   % in-phase part P (or R)   I = imag(cross);    % quadrature part Q (or I)   R = R(2:T/2+1).^2;  % only the first half is used % R = R.*2; I = I(2:T/2+1).^2; % only the first half is used % I = I.*2;

xhat = abs(xhat); xhat = xhat(2:T/2+1).^2; % only the first half is used % xhat = xhat.*2;

yhat = abs(yhat); yhat = yhat(2:T/2+1).^2; % only the first half is used % yhat = yhat.*2;

scs = (R+I)./(xhat.*yhat); % equal to 1!!!

%*******************************************************  % HOWEVER, if all the values are used:                 * % R = R.*2; I = I.*2; xhat = xhat.*2; yhat = yhat.*2; * % scs = (R+I)./(xhat.*yhat); % EQUAL to 0.5, NOT 1!!! *  %*******************************************************

plot(f,scs,'-k','LineWidth',1.0); xlabel('frequency (cycle/month)','FontSize',12); ylabel('coherency','FontSize',12); grid on; ylim([0 2]); set(gca,'FontSize',12); title('Squared Coherency Spectrum','FontSize',12); print('-djpeg','ch_4.5') elseif flag==6 %__________________________________________________________________________  % 4.6 squared coherency spectrum after rebinning f = [1:T/2]'/T;

xhat = fft(ts1)/T; yhat = fft(ts2)/T; cross = xhat .* conj(yhat); R = real(cross);   % in-phase part P (or R)   I = imag(cross);    % quadrature part Q (or I)   R = R(2:T/2+1).^2; I = I(2:T/2+1).^2;

xhat = abs(xhat); xhat = xhat(2:T/2+1).^2; yhat = abs(yhat); yhat = yhat(2:T/2+1).^2;

j = 0; for i=1:3:117 j = j+1; fn(j) = (f(i)+f(i+1)+f(i+2))/3; Rn(j) = (R(i)+R(i+1)+R(i+2))/3; In(j) = (I(i)+I(i+1)+I(i+2))/3; xhatn(j) = (xhat(i)+xhat(i+1)+xhat(i+2))/3; yhatn(j) = (yhat(i)+yhat(i+1)+yhat(i+2))/3; end scs = (Rn+In)./(xhatn.*yhatn); % not equal to 1 anymore

plot(fn,scs,'-k','LineWidth',1.0); xlabel('frequency (cycle/month)','FontSize',12); ylabel('coherency','FontSize',12); grid on; axis([fn(1) fn(end) 0 3]); set(gca,'FontSize',12); title('Squared Coherency Spectrum after rebinning','FontSize',12); print('-djpeg','ch_4.6') elseif flag==7 %__________________________________________________________________________  % 4.7 phase spectrum f = [-(T-1)/2:T/2]/T;

xhat = fft(ts1)/T; yhat = fft(ts2)/T; cross = xhat .* conj(yhat); R = real(cross);   % in-phase part P (or R)   I = imag(cross);    % quadrature part Q (or I)   R = cat(1,flipud(R(122:end)),R(1:121)); R = cat(1,flipud(I(122:end)),I(1:121));

pspct = atand(I./R);

plot(f,pspct,'x-','color',[0.3 0.1 0.5],...      'MarkerSize',5,'MarkerEdgeColor',[0.1 0.5 0.1],'LineWidth',1.0); xlabel('frequency (cycle/month)','FontSize',12); ylabel('phase','FontSize',12); grid on; ylim([-90 90]); set(gca,'FontSize',12); daspect([1 1000 1]); title('Phase Spectrum','FontSize',12); print('-djpeg','ch_4.7') end code

**Question 5**
code format="matlab" % Specify 'flag' to do part 5.2 or 5.3 clear;clc;close all; flag = 3;   % 2: part 5.2; 3: part 5.3

% grab data: sst, slp sst = nc_varget('data.nc','sst'); slp = nc_varget('data.nc','slp');

T = 240; f = [-119:120]/T;  % frequency goes from -119/240 to 120/240 L = 144; k = [-71:72]/L;    % wavenumber goes from -71/114 to 72/114

if flag==2 %__________________________________________________________________________  % 5.2 2D fft p = fft2(sst)/(T*L);     % the same as fft(fft(sst).').'; p = abs(p); p = fftshift(p); % p(121,73)=0; contourf(k,f,log10(p),1000,'lines','none') % log(p)!!! % caxis([0 0.1]); % caxis([-2.5 0.5]) grid on; set(gca,'FontSize',12); xlabel('wavenumber (2.5^o ^-^1)','FontSize',12); ylabel('frequency (month^-^1)','FontSize',12); title('2D Power Spectrum','FontSize',12); set(gca,'FontSize',12); print('-djpeg','ch_5.2.log') else %__________________________________________________________________________  % 5.3 filtered data image p = fft2(sst); p = fftshift(p); f(f>1/108 & f<1/24)=0; f(f<-1/108 & f>-1/24)=0; f(f~=0)=1; % band pass I: removes a set of low freq. % f(f>2/12 & f<6/12)=0; f(f<-2/12 & f>-6/12)=0; f(f~=0)=1; % band pass II: removes a set of high freq. for i=1:144 p(:,i)=p(:,i).*f'; end p = fftshift(p); sst_f = ifft2(p);

subplot(1,2,2); contourf(abs(sst_f)); daspect([1 1.5 1]); caxis([20 30]); colorbar('location','SouthOutside'); xlabel('Longitude','FontSize',14,'FontWeight','bold'); ylabel('Time/month','FontSize',14,'FontWeight','bold'); % title('SST after band pass II','FontSize',14,'FontWeight','bold'); title('SST after band pass I','FontSize',14,'FontWeight','bold'); set(gca,'XTick',[1 24 48 72 96 120 144]); set(gca, 'XTickLabel',{' 0';' 60';'120';'180';'240';'300';'360'},...       'FontSize',14,'FontWeight','bold');

subplot(1,2,1); contourf(sst; daspect([1 1.5 1]); caxis([20 30]);  xlabel('Longitude','FontSize',14,'FontWeight','bold');   ylabel('Time/month','FontSize',14,'FontWeight','bold');   title('SST original','FontSize',14,'FontWeight','bold');   set(gca,'XTick',[1 24 48 72 96 120 144]);   set(gca, 'XTickLabel',{'  0';' 60';'120';'180';'240';'300';'360'},... 'FontSize',14,'FontWeight','bold');  print('-djpeg','ch_5.3') end

code