【物理应用】基于matlab波数谱计算【含Matlab源码 508期】

一、获取代码方式

获取代码方式1：
完整代码已上传我的资源：【物理应用】基于matlab波数谱计算【含Matlab源码 508期】

获取代码方式2：
通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。

备注：
订阅紫极神光博客付费专栏，可免费获得1份代码（有效期为订阅日起，三天内有效）；

二、部分源代码

clear
figure
load tide36part16RSPO
%load tide36all

load all_atm_pres.mat
a=mean(atmp');
pa=(a-mean(a))'/100;
t=[166.5:0.5:901.5]';
%t=t(2:2:1472);
%pa=interp1(t,pa,onetime1n);
tidem1d5B=interp1(time36part,tide36part,t);
tidem=tidem1d5B;

tidem1d5pa=tidem+pa(2:1:1472);


load uriba5
ba=interp1(onetime,BA5,t);

aa=tidem1d5pa+ba;
st=isnan(aa);
ba(st)=[];
tidem1d5pa(st)=[];
t(st)=[];

in=ba;order=4;filtT=[70];sampT=1/2;
%hold off
%plot(ba);hold on;plot(ba-ba70','r')

filtT=[70 2];
filtT=[70];

%ba70n=butterfilthigh(in, order,70,sampT)*100;
%ba70=butterfiltbandpass(in, order,filtT,sampT)*100;
ba70=butterfilthigh(in, order,filtT,sampT)*100;

%plot(ba70,'g')
%plot(ba);hold on;plot(ba-ba70','r')
in=tidem1d5pa;
%tide70=butterfiltbandpass(in, order,filtT,sampT)*100;
tide70=butterfilthigh(in, order,filtT,sampT)*100;

nFFT=128*2;
Fs=2;
WinLength=nFFT;
nOverlap=WinLength/2;
Window=hanning(WinLength); %Window=(WinLength);
%Window=WinLength;
Detrend='linear';
   P=0.90;
 subplot(2,2,3);
    hold off

ecp2=ba70;  
[psd_bp1,pconf,f]=psd(ecp2,nFFT,Fs,Window,nOverlap,P);
semilogx(f,psd_bp1.*f,'k');
xlim([0.01 1])

xlabel('cycles/day');hold on
ylabel('cm^2');hold on

 ecp2=tide70;  
[psd_bp1,pconf,f]=psd(ecp2,nFFT,Fs,Window,nOverlap,P);
semilogx(f,psd_bp1.*f,'r','color',[0.5 0.5 0.5]);
xlabel('cycles/day');
title('Variance preserving PSD')

%NW=fix(736/64*2); x=tide70;y=tide70;qbias=[];confn=[];qplot=1; dt=1;
%[s, c, ph, ci, phi] = cmtm(x,y,dt,NW,qbias,confn,qplot);ci(1)

%h=legend ('Tide station','Bottom pressure');
%h1=findobj(h,'type','text')
%set(h1(1),'color',[0.5 0.5 0.5])
%set(h1(2),'color',[0 0 0])


subplot(4,2,6)
[ Pxx, Pyy, Pxy, coh, pha, freq ] = func_coherence((ba70),(tide70), nFFT, Fs, nFFT, nOverlap );
 f=freq;
semilogx(f,coh,'k');hold on
%plot(f,ones(65,1)*0.3730,'k')
plot(f,ones(129,1)*0.3730,'k')
function [out]=butterfilt(in,order,filtT,sampT);
%
% function [out]=butterfilt(in,order,filtT,sampT);
%
% This function is designed to do lowpass filtering of any vector of 
% data.  It uses a  Butterworth filter and does the 
% filtering twice, once forward and once in the reverse direction.  
%
% INPUTS:
%        in    -- the data vector
%	 filtT -- the smoothing period; 
%	 sampT -- the sample period
%        order -- of the butterworth filter
%
%  NOTE:  filtT and sampT must be entered in the same units!!!!!
%

[mm,nn]=size(in);
if nn>mm & mm==1
 in=in';
 [mm,nn]=size(in);
elseif mm~=1 & nn~=1
 error('Input must be a vector and not an array!!!')
end

%
% First remove a ramp so that the first and last points are zero.  
% This minimizes the filter transients at the ends of the record.  
%
P=polyfit([1 mm],[in(1) in(mm)],1);
dumx=[1:mm]';
lineartrend=polyval(P,dumx)*0;
in=in-lineartrend;

%
% Second, create the Butterworth filter.  
%
Wn=(2./filtT)*sampT;
[b,a]=butter(order,Wn,'high');

%
% Third, filter the data both forward and reverse.  
%
inan=find(isnan(in));
if ~isempty(inan)
 error('There are NaNs in the data set!!! Remove them and try again!')
end
out=filtfilt(b,a,in);
%

三、运行结果

四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 门云阁.MATLAB物理计算与可视化[M].清华大学出版社，2013.

文章来源: qq912100926.blog.csdn.net，作者：海神之光，版权归原作者所有，如需转载，请联系作者。

原文链接：qq912100926.blog.csdn.net/article/details/114707723