FKT变换程序手册

F-K变换程序手册编译环境：MatlabR2013a64-bit(8.1.0.604)编译时间：2017.6.20源程序列表序号文件名称说明1fktran.mF-K变换2fftrl.m傅里叶变换3demo.m主函数4R.mat调用文件5t.mat调用文件6x.mat调用文件1、fktran.mF-K变换程序说明：F-K变换函数原型：function[spec,f,kx]=fktran(seis,t,x,ntpad,nxpad,percent,ishift)参数列表：序号参数名说明1seis输入二维地震记录2nsamp,ntr输入信号维度3nxpad将位置序列补充为2的整数次幂4ntpad将时间序列补充为2的整数次幂5spec频率波数谱6kx波数序列7f频率序列8t输入时间序列9x输入位置序列附程序代码：function[spec,f,kx]=fktran(seis,t,x,ntpad,nxpad,percent,ishift)%FKTRAN:Forwardfktransform%[spec,f,kx]=fktran(seis,t,x,ntpad,nxpad,percent,ishift)%%FKTRANusesMatlab'sbuiltinffttoperforma2-df-ktransform%onarealvalued(presumablyseismictime-space)matrix.Onlythe%positivef'sarecalculatedwhileallkxsare.Matlab's'fft'isused%oneachcolumnforthet-ftransformwhile'ifft'isusedoneachrow%forthex-kxtransform.Theinversetransformisperformedbyifktran.%%seis...input2-dseismicmatrix.Onetracepercolumn.%t...vectoroftimecoordinatesforseis.length(t)mustbethe%sameasnumberofrowsinseis.%x...vectorofspacecoordinatesforseis.length(x)mustbethesame%asthenumberofcolumnsinseis.%ntpad...padseiswithzerofilledrowsuntilitisthissize.%*******default=nextpowerof2******%nxpad...padseiswithzerofilledcolumnsuntilitisthissize.%*******default=nextpowerof2******%percent...applyaraisedcosinetapertobothtandxpriortozeropad.%lengthoftaperistheispercentageofthelengthofthexandtaxes%*******default=0*********%ishift...if1,thenthekaxisofthetransformisunwrappedtoput%kx=0inthemiddle.%*******default=1*******%spec...complexvaluedf-ktransformofseis%f...vectoroffrequencycoordinatesfortherowsofspec%kx...vectorofwavenumbercoordinatesforthecolumnsofspec%[nsamp,ntr]=size(seis);if(length(t)~=nsamp)error('Timecoordinatevectorisincorrect');endif(length(x)~=ntr)error('Spacecoordinatevectorisincorrect');endif(nargin<7)ishift=1;endif(nargin<6)percent=0.;endif(nargin<5)nxpad=2^nextpow2(x);endif(nargin<4)ntpad=2^nextpow2(t);end%usefftrlforthet-ftransform[specfx,f]=fftrl(seis,t,percent,ntpad);clearseis;%oktaperandpadinxif(percent>0)mw=mwindow(ntr,percent);mw=mw(ones(1,ntr),:);specfx=specfx.*mw;clearmw;endif(ntr<nxpad)ntr=nxpad;%thiscausesiffttoapplythexpadend%fftonrowsspec=ifft(specfx.',ntr).';%computekxkxnyq=1./(2.*(x(2)-x(1)));dkx=2.*kxnyq/ntr;kx=[0:dkx:kxnyq-dkx-kxnyq:dkx:-dkx];if(ishift==1)[kx,ikx]=sort(kx);spec=spec(:,ikx);end2、fftrl.m傅里叶变换程序说明：傅里叶变换函数原型：function[spec,f]=fftrl(s,t,percent,n)参数列表：序号参数名说明1s输入信号2t时间序列3spec频谱4f频率附程序代码：function[spec,f]=fftrl(s,t,percent,n)%FFTRL:forwardFouriertransformforreal-valuedsignals.%%[spec,f]=fftrl(s,t,percent,n)%[spec,f]=fftrl(s,t,percent)%[spec,f]=fftrl(s,t)%%Forwardfouriertransformofareal-valuedsignal.Relativeto%MATLAB'sfftitreturnsonlythepositivefrequenciesinanarray%roughlyhalfthesize.Ifthereareninputtimesamples,thenthere%willbefloor(n/2+1)frequencysamples.Thismeansthatifnisaneven%number,thenatimeseriesoflengthnandoneoflengthn+1willproduce%frequencyspectraofthesamelength.However,onlythefirstcasewill%haveasampleatNyquist.%Iftheinputtraceisavector,thenthereturn%issimplythetransformofthatvector.Ifamatrix,theneach%columnofthematrixistransformed.Theinverseisaccomplishedby%IFFTRL.%%s=inputsignal(traceorgather,i.e.vectorormatrix)%t=inputtimecoordinatevector%percent=specifiesthelengthofaraisedcosinetapertobe%appliedtos(bothends)priortoanypadding.%Taperisapercentofthelengthofs.Taperis%appliedusingMWINDOW.%**********Default=0%***********%n=lengthtowhichtheinputtraceistobepaddedwithzeros.%**********Defaultistheinputlength(nopad)***********%%spec=outputspectrum%f=outputfrequencysamplevectorif(nargin<4)n=length(t);endif(nargin<3)percent=0.0;end%determinenumberoftracesinensemble[l,m]=size(s);ntraces=1;itr=0;%transposeflagif(l==1)nsamps=m;itr=1;s=s(:);%switchtocolumnvectorselseif(m==1)nsamps=l;elsensamps=l;ntraces=m;endif(nsamps~=length(t))t=t(1)+(t(2)-t(1))*(0:nsamps-1);if(nargin<4)n=length(t);end%error('timevectorandtracematrixdon''tmatchinlength');end%applythetaperif(percent>0)mw=mwindow(nsamps,percent);mw=mw(:,ones(1,ntraces));s=s.*mw;clearmwend%padsifneededif(nsamps<n),s=[s;zeros(n-nsamps,ntraces)];nsamps=n;end%transformthearrayspec=fft(s,nsamps);spec=spec(1:floor(n/2+1),:);%saveonlythepositivefrequenciesclears;%buildthefrequencyvectorfnyq=1./(2*(t(2)-t(1)));nf=size(spec,1);df=2*fnyq/n;%f=linspace(0.,fnyq,nf)';f=df*(0:nf-1)';if(itr)f=f';spec=spec.';end3、demo.m主函数程序说明：主函数