slope_adapt_filter.m

函数实现对InSAR复相位信号的基于局部频率估计的“斜坡自适应滤波” 此滤波方法特别适用于InSAR相位图滤波

function phasef=slope_adapt_filter(phase,xsize,ysize)
%基于局部频率估计的斜坡自适应滤波（利用FFT+CZT法提高了局部频率估计的精度）
windowsizei=13;%估计窗方位向尺寸
windowsizej=13;
pwsizei=13;%滤波窗方位向尺寸
pwsizej=13;
wini=windowsizei;
wi=(wini-1)/2;
winj=windowsizej;
wj=(winj-1)/2;
pwhi=(pwsizei-1)/2;
pwhj=(pwsizej-1)/2;
phaseb=zeros(xsize+wini-1,ysize+winj-1);
dimb=size(phaseb);
Nr=dimb(1);
Na=dimb(2);
phaseb(wi+1:Nr-wi,wj+1:Na-wj)=phase;
phaseb(1:wi,1:Na)=flipud(phaseb(wi+1:2*wi,1:Na));
phaseb(1:Nr,1:wj)=fliplr(phaseb(1:Nr,wj+1:2*wj));
phaseb(Nr-wi+1:Nr,1:Na)=flipud(phaseb(Nr-2*wi+1:Nr-wi,1:Na));
phaseb(1:Nr,Na-wj+1:Na)=fliplr(phaseb(1:Nr,Na-2*wj+1:Na-wj));
phaseb=exp(j*phaseb);%将相位转化为复数
phi=zeros(Nr,Na);
phasef=zeros(Nr,Na);
nn=32;%FFT2 的点数
mm=32;%CZT变换的点数mm
mn=mm*nn;%频谱总分辨率（mn×mn）
phasej=zeros(windowsizei,windowsizej);
tempi=0:pwsizei-1;
tempj=(0:pwsizej-1)';
w0=exp(-j*2*pi/mn);
for m=wi+1:Nr-wi
  for n=wj+1:Na-wj
     phasej(1:wini,1:winj)=phaseb(m-wi:m+wi,n-wj:n+wj);
     phasefj=fft2(phasej,nn,nn);%窗口内二维傅立叶
     phasefj=fftshift(phasefj);
     [C,index_x]=max(max(abs(phasefj)));%index得到的是最大值对应的列序数
     [D,index_y]=max(abs(phasefj(:,index_x)));%index得到的是最大值对应的行序数
     p=double(index_x-nn/2-2)/nn;
     q=double(index_y-nn/2-2)/nn;
     ex=exp(j*2*pi*p);
     ey=exp(j*2*pi*q);
     %下面针对极值点进行CZT变换主瓣频谱精细化
     phasef_czt1=czt(phasej,3*mm,w0,ey);%方位向czt
     phasef_czt2=czt(phasef_czt1.',3*mm,w0,ex);%距离向czt
     phasef_czt2=phasef_czt2.';
     [C,index_x]=max(max(abs(phasef_czt2)));
     [D,index_y]=max(abs(phasef_czt2(:,index_x)));
     pp=double(index_x-1)/mn;%细化后偏移量
     qq=double(index_y-1)/mn;
     ppp=p+pp;%在该分辨率整个频谱中的位置
     qqq=q+qq;
     ex=exp(j*2*pi*ppp);%距离向局部频率（指数的）
     ey=exp(j*2*pi*qqq);%方位向局部频率（指数的）
     phi(m-pwhi:m+pwhi,n-pwhj:n+pwhj)=ey.^tempj*(ex.^tempi);%窗口线性相位部分
     phish=phaseb(m-pwhi:m+pwhi,n-pwhj:n+pwhj).*conj(phi(m-pwhi:m+pwhi,n-pwhj:n+pwhj));
     phasef(m,n)=mean(mean(phish))*phi(m,n);%均值滤波,得到经局部频率线性平滑的相位
  end
  waitbar(m/Nr)
end
phasef=angle(phasef);
phasef=phasef(wi+1:Nr-wi,wj+1:Na-wj);%取有效范围