qiyizhi_潮汐奇异值_样条插值_潮汐_

clc; clearvars -except XM; %导入文件变量的过程是用matlab的变量功能实现的，导入的变量作为数值矩阵XM存放在工作区内。 zz=transpose(XM); zzz=zz(:); a=zzz(:)';%将所有数值矩阵XM所有的值转化为行向量，方便作为纵坐标进行绘图，从而对奇异值进行初筛。 b=1:1:744; %因为奇异值为函数图像的极大值或极小值，首先第一步找出图像的极值点。 %计算原序列的平均值和标准差。 ave1=mean(a); s1=sum(((a(1,:)-mean(a)).^2)/(length(a)-1)).^0.5; disp(['原序列的平均值为',num2str(ave1)]); disp(['原序列的标准差为',num2str(s1)]); for i=2:1:743 %找出极大值点的横坐标并赋值给max if (a(i)>a(i-1) && a(i)>a(i+1)) max(i)=i; end end max(max==0)=[]; for i=2:1:743 %找出极小值的横坐标并赋值给min。 if (a(i)<a(i-1) && a(i)<a(i+1)) min(1,i)=i; end end min(min==0)=[]; %在找完极值点后，就要对奇异值进行筛选，通过观察函数图像可以得出一个比较明显的结论， %即所有的奇异值点的左侧或右侧的第一个数据点也必为极值点。这是因为奇异值点为数据不准确 %的点，潮汐数据在涨落时间内水位变化是单向变化，那么如果在某一时间段内出现反向变化的点， %那就是潮汐数据的奇异值。即奇异值两侧至少有一个数据点为极值点。 max1=[]; for c=1:1:64 if (a(max(c)-1)-a(max(c)-2)<0 || a(max(c)+1)-a(max(c)+2)<0) max1=[max1,max(c)]; end end min1=[]; for d=1:1:61 if (a(min(d)-1)-a(min(d)-2)>0 || a(min(d)+1)-a(min(d)+2)>0) min1=[min1,min(d)]; end end %max1和min1就是初筛后的有可能是奇异点的数据，max对应的是极大值点，min对应的是极小值点。 %接下来就要进行二次筛选。 %对于极大值点可分为两种情况进行讨论，如果左侧为极值点，即说明曲线是下降趋势，那么就要满足 %左侧数据大于右侧数据（因为左侧是极小值点），才能保证去掉奇异值后该段曲线仍然成下降趋势；反之，若右侧为极值点，则为上升趋势，则要满足 %左侧数据小于右侧数据。同时，因为奇异值两侧的点需要满足一定的线性趋势，因此，需要加上限定条件。 %通过观察潮汐数据，发现相邻的数据点只差的绝对值大致小于300（可由算法验证，较为简单，所以程序里没有体现）， %可由此作为限定条件进行严格筛选。 max2=[];max3=[]; for n=1:4 if (a(max1(n)-1)-a(max1(n)-2)<0 && a(max1(n)+1)-a(max1(n)-1)<-300) max2=[max2,max1(n)]; end if (a(max1(n)+1)-a(max1(n)+2)<0 && a(max1(n)-1)-a(max1(n)+1)<-300) max3=[max3,max1(n)]; end end maxm=[max2,max3]; %筛选出第508和604个数据点为奇异值点。 min2=[];min3=[]; %接着筛选极小值点内的奇异值，思路同上。 for n=1:4 if (a(min1(n)+1)-a(min1(n)+2)>0 && a(min1(n)-1)-a(min1(n)+1)>300) min2=[min2,min1(n)]; end if (a(min1(n)-1)-a(min1(n)-2)>0 && a(min1(n)-1)-a(min1(n)+1)<-300) min3=[min3,min1(n)]; end end minm=[min2,min3]; %筛选出第229和386个数据点为奇异值点。 qiyi=[maxm,minm]; disp('以下几个数对应的数据点为奇异值'); disp(qiyi); %将奇异值点的横坐标赋值给矩阵qiyi，以方便后面的插值计算，同时在屏幕上输出。 %我一共尝试了三种插值方法，发现除了二次插值的准确度较差以外，线性插值和样条插值的准确度 %都比较高，故最后保留了样条插值的代码，并将插值后的数据点存放在一个新的txt文件里。 %若要尝试其他的插值方法，可直接调用下面的函数。 new=XM; for n=1:4 x=[qiyi(n)-1 qiyi(n)+1]; y=[a(qiyi(n)-1) a(qiyi(n)+1)]; xx=qiyi(n)-1:0.1:qiyi(n)+1; yy=spline(x,y,xx); a(qiyi(n))=yy(11); new(qiyi(n))=yy(11); end file=fopen('f:\data.txt','w'); fprintf(file,'%.0f ',new); sta=fclose(file); plot(b,a); %将插值替代后的数据作为纵坐标，进行绘图。 %计算插值替代数据后的均值和标准差，并在屏幕上输出对比 ave2=mean(a); s2=sum(((a(1,:)-mean(a)).^2)/(length(a)-1)).^0.5; disp(['新序列的平均值为',num2str(ave2)]); disp(['新序列的标准差为',num2str(s2)]); function [l]=lagrange(x,x0,x1,y0,y1) %线性插值函数 l0=(x-x1)/(x0-x1); l1=(x-x0)/(x1-x0); l=l0*y0+l1*y1; end function [r]=erci(x,x0,x1,x2,y0,y1,y2) %二次插值函数 l0=(x-x1)*(x-x2)/(x0-x1)*(x0-x2); l1=(x-x0)*(x-x2)/(x1-x0)*(x1-x2); l2=(x-x0)*(x-x1)/(x2-x0)*(x2-x1); r=y0*l0+y1*l1+y2*l2; end