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

标签：cycle   NPU   mil   filter   trend   %s   smooth   set   center   

一、简介

基于matlab波数谱计算

二、源代码

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

三、运行结果

四、备注

版本：2014a
完整代码或代写加1564658423

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

标签：cycle   NPU   mil   filter   trend   %s   smooth   set   center   

原文地址：https://www.cnblogs.com/homeofmatlab/p/14942696.html