clear,close all;
clc;
load fk2013meteor_equinoxMon9and10.mat;
%% Select the height
% height_bin=Height_Year_equinox(:,1,1);
height_bin=70:2:110;
aa=88:2:98;
bb=ismember(height_bin,aa);
%% get the Meridional Wind (MW) and Zonal Wind (ZW)
cc=MeridianVel_Year_equinox(bb,:,:); % MW
dd=ZonalVel_Year_equinox(bb,:,:); % ZW
ee_Z90=dd(2,:,:);
ff_Z90=squeeze(ee_Z90);
gg_Z90=ff_Z90(:); % MW90
%%
sst = gg_Z90; % 1464 X 1
%%
%% interpolate with Spline method
% calculate the max, min and mean
sst_min=min(sst);
sst_max=max(sst);
sst_nanmean=nanmean(sst);
% interpolate with Spline
X=find(isnan(sst)==0);
x=find(isnan(sst)==1);
z=interp1(X,sst(X),x,'spline');
n=length(x);
for iii=1:n
sst(x(iii))=z(iii);
end
% remove singularities
sst(sst>sst_max | sst<sst_min)=[sst_nanmean];
sst_interp=sst;
%% band-pass filtered, mothod 1-zhaoguangxin
sst=sst_interp;
sst_len=length(sst);
Wn=[2/7/24 2/6/24]; % 1-16 day band-pass filtered
[B,A]=butter(4,Wn);
VX=[sst;sst;sst];
VX2=filtfilt(B,A,VX);
VX_output=VX2(sst_len+1:2*sst_len);
sst_bpf1=VX_output;
%% band-pass filtered, mothod 2-zhuzhengping
x1=sst_interp;
fs=24;
fp1=(6.5/(6.5+1))*1/6.5;
fp2=(6.5/(6.5-1))*1/6.5;
fs1=(6.5/(6.5+1+1))*1/6.5;
fs2=(6.5/(6.5-1-1))*1/6.5;
Wp1=[fp1*2*pi/fs,fp2*2*pi/fs]./pi;
Ws1=[fs1*2*pi/fs,fs2*2*pi/fs]./pi;
Ap1=3;
As1=40;
% [N1,Wn1]=buttord(Wp1,Ws1,Ap1,As1);
% [b1,a1]=fir1(N1,Wn1);%,'bandpass');
% y1=filter(b1,a1,x1);
% x1=[x1;x1;x1];
ff = design(fdesign.bandpass(Ws1(1),Wp1(1),Wp1(2),Ws1(2),As1,Ap1,As1,fs),'butter');
y1=filter(ff,x1);
% y2=y1(sst_len+1:2*sst_len);
subplot(3,1,1);
plot(sst_interp,'r*-');
title('Original data');
% hold on;
subplot(3,1,2);
plot(sst_bpf1,'go-');
title('Method 1:zhaoguangxin method');
subplot(3,1,3);
plot(y1,'bd-');
title('Method 2:zhuzhengping method');
% sst_len=length(sst);
% fs=1; % 1 point/hour
%% lomb algorithm
% timeday02=0:length(sst)-1;
% timeday02=timeday02';
% % ffinterpolate = sst; % [2908 X 1], data with the nans,
% % datalength=size(ffinterpolate);
%
% frequency_Day=(0:1/500:1/2)./24; % Period is from 500-day to 2-day
% period_Day=(1./frequency_Day(1:25:end))./24; % period in day
% period_Day_str=sprintf('%4.2f %20.2f %18.2f %18.2f %20.2f %20.2f %20.2f %20.2f %20.2f %20.2f %20.2f',period_Day(1:end));
%
% %%%
% [Px, Probx] = lomb_zzp_fun(timeday02, sst, frequency_Day); %%Px is amplitude
% % [Py, Proby] = lomb_zzp_fun(Ty, Vy, freq);
%
% %%
% [pks,locs] = findpeaks(Px);
%
% pks
% ooo=frequency_Day(locs)
%
% %%%%%%%%%%%%%%%%%%%%%%
% hf=figure(1);
%
% subplot(2,1,1)
% plot(frequency_Day,Px,'r.-')
% set(gca,'xlim',[min(frequency_Day),max(frequency_Day)]);
% set(gca,'xtick',frequency_Day(1:25:end),'xticklabel',frequency_Day(1:25:end));
%
% xlabel('Frequency (cycles/day)');
% ylabel('Amplitude (m/s)');
% title('spectrum of zonal wind at 90 km');
% aaaa=get(gca,'ylim');
% text(0,aaaa(2)-1.5,period_Day_str);
% subplot(2,1,2)
% Probx_positive=1-Probx;
% plot(frequency_Day,Probx_positive*100,'g.-');
% set(gca,'xlim',[min(frequency_Day),max(frequency_Day)],'ylim',[0 110]);
%% print to the file
% set(hf, 'Renderer', 'painters');
% set(hf,'PaperPositionMode','auto');
% figname=[pwd '\' 'figures\fig97_sy2013_equM910_lomb_Z90.tif'];
% set(hf, 'PaperUnits', 'centimeters');
% % set(gcf, 'PaperPosition', [0 0 6 12]);
% set(hf, 'PaperPosition', [0 0 29 18]);
% % set(gcf, 'PaperPosition', [0 0 38 24]);
% print(hf, '-dtiff','-r300',figname);
