频率切片小波变换

function[Y1]=GetInvFSWT(N,a,fp)%%fp:Observationfrequency[Tn,np]=size(a);%%Tn:TimeResamplePointNumber,npisthelengthoffp;%%%nl=length(fp);%%Observationfrequencylengths1=zeros(N,1);forp=1:np;sum=0;fork=1:Tn;%%ifmod(k-1,2)==0sum=sum+conj(a(k,p));elsesum=sum-conj(a(k,p));endsum=sum+conj(a(k,p));endsum=sum/Tn;iffp(p)>0s1(fp(p)+1)=sum;s1(N-fp(p)+1)=conj(sum);endends1(1)=0;Y1=fft(conj(s1),N)/N;%Y1=ifft(s1,N);Y1=real(Y1);function[a]=GetFSWT(s,Fs,fp,kapa,Tn)%%Tn:TimeResamplePointNumber%%fp:ObservedfrequencyrangeN=length(s);%%%OriginalPointNumberY=fft(s,N);%%FFTY(1)=0;%%cuttheconstantpartnl=length(fp);%%Observationfrequencylengtha=zeros(Tn,nl);s1=zeros(Tn,1);TNN=fix(Tn/2);k0=0;forp=1:nl;s1=s1*0;k0=0;fork=1:Tn;s1(k)=0;iffp(p)==0continue;end;m=fp(p)+k-1-TNN;ifm>=0&&abs(k-1-TNN)<.5*TNN&&m<Nct=kapa*(k-1-TNN)/fp(p);% s1(k)=Y(m+1)/(1+ct*ct);s1(k)=Y(m+1)*exp(-ct*ct/2); %%%HereyoucanuseyourFrequencyslicefunction.end;ends1(1)=0;s1=conj(s1);Y1=fft(s1,Tn);fork=1:Tn;ifmod(k-1,2)==1Y1(k)=-Y1(k);end;enda(:,p)=Y1;end测试：t=(0：1000)/1000；s=sin(2*pi*10*t)+0.7*sin(2*pi*28*t+pi/3);Fs=1000;%%采样率N=length(s);s=s-sum(s)/N;%%PleasecuttheDCpartinthesignal%%s0=s;s=s/max(s)+sin(rand(N,1)*2*pi)*0.5;%%s=s.*(1+sin(rand(N,1)*2*pi));f1=0;%%youcanchangeitf2=200e6;%%youcanchangeit%%%[f1,f2]isyourobservedfrequencyrange,youcanchangeanyscopek1=fix(f1*N/Fs-0.5);k2=fix(f2*N/Fs-0.5);df=1;%%isaobservedfrequencysteplengthif(k2>N/2+1)k2=N/2+1;endfp=fix(k1:df:k2);%%fpiftheobservedfrequencyindiscreteformnl=length(fp);kapa=sqrt(2)/2/0.025;%%kapaisthetime-frequencyresolutionfactorTn=512;%%smaplenumerintimedomain,youcanchageit[A]=GetFSWT(s,Fs,fp,kapa,Tn);%获取指定区域时频谱%A1=zeros(size(A));%A1(80:160,400:650)=A(80:160,400:650);%A=A1;A=A;ss=GetInvFSWT(N,A,fp);%recoverthesignalfromthebandoffp;B=sqrt(A.*conj(A));mx=max(max(B));B=fix(B*128/mx);t=(0:Tn-1)*N/Fs/Tn;figure(1)subplot(221)plot((0:N-1)/Fs,s,'b');holdonplot((0:N-1)/Fs,ss,'r');subplot(222)Y=fft(s,N);z=Y.*conj(Y);z=sqrt(z);z(1)=0;z(2)=0;K1=fix(f2*N/Fs);fp1=[0:K1-1]/N*Fs;plot(fp1,z(1:K1));subplot(223)

%image(f