filter_bp_by_Zhaogz_Zhuzp.rar_bandpass_filter_growewe_orange7rz_

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);