风生近惯性振荡，用surface drifter数据滤波，得到近惯性流速 matlab

clear,clc,close all
%%
cd D:\y_data\GDP\gps
filelist = dir(‘D:\y_data\GDP\gps*.nc’);
k = length(filelist);
omega=7.292e-5;
inertial_current=[];
x=1;
for i = 1:k
filename=filelist(i).name;
lat = ncread(filename,‘latitude’);
lon = ncread(filename,‘longitude’);
u = ncread(filename,‘ve’);
v = ncread(filename,‘vn’);
time = ncread(filename,‘time’);
date = datestr(time/86400+datenum(1970,1,1),‘yyyymmdd’);
if (min(lat)>=10 & max(lat)<=70) & (min(lon)>=-90 & max(lon)<=30)
f = abs(sw_f(lat));
dt=3600;fs=24; %采样频率
if length(time)>300
for j =1:300:length(f)
try
f1=f(j:j+299);
u1=u(j:j+299);
v1=v(j:j+299);
time1=time(j:j+299);
lon1=lon(j:j+299);
lat1=lat(j:j+299);
catch
f1=f(end-299:end);
u1=u(end-299:end);
v1=v(end-299:end);
time1=time(end-299:end);
lon1=lon(end-299:end);
lat1=lat(end-299:end);
end
diff_time=diff(time1);
if all(diff_time==3600)
fc1=mean(0.75f)360024/(2pi);
fc2=mean(1.25f)360024/(2pi);
% if fc2>fs/2
% fc2=1.388e-4;
% end
if fc2 < fs/2
wn=[fc1 fc2]/(fs/2);
[b a]=butter(4,wn);
udd=filtfilt(b,a,double(u1));
vdd=filtfilt(b,a,double(v1));
% u1=double(u1);
% v1=double(v1);
% udd=bandpass2(u1’,[fc1 fc2],fs,6);
% vdd=bandpass2(v1’,[fc1 fc2],fs,6);
inertial_current(x).time=time1;
inertial_current(x).lat=lat1;
inertial_current(x).lon=lon1;
inertial_current(x).udd=udd;
inertial_current(x).vdd=vdd;
inertial_current(x).u1=u1;
inertial_current(x).v1=v1;
x=x+1;
end
end
end
end
end
end

每300个小时做一个滤波
频率在0.75f-1.25f