%%%?二維傅里葉變換相位提取

clear?all

N=512;??????%?抽樣點(diǎn)數(shù)
L=1;????????%?x，y方向長(zhǎng)度
a=1;????????%?背景強(qiáng)度
b=0.8;??????%?調(diào)制度
r0=1/16;
f_x=r0*N;????%?x方向載頻
f_y=r0*N;????%?y方向載頻
delta=L/N;???%?空域采樣間隔
delta_f=1/（N*delta）;?%?頻域采樣間隔

[x?y]=meshgrid（（-N/2:N/2-1）*delta）;
[fx?fy]=meshgrid（（-N/2:N/2-1）*delta_f）;????

[XY]=meshgrid（1:N）;

W=0.65e-6.*（x.^2+y.^2）.^2;
PV=max（max（W））-min（min（W））;
RMS=sqrt（sum（sum（W.^2）））/N;
lamda=532e-9;
k=2*pi/lamda;
phi_real=k.*W;
i=a+b.*cos（2*pi*f_x.*x+2*pi*f_y.*y+phi_real）;
bg=sum（sum（i））/（N*N）;
??
?i1=uint8（255*i）;
?figureimagesc（i1）;
?colormap?gray
?axis?square
?
flag_1=（（sqrt（（X-N/2）.^2+（Y-N/2）.^2））>N/2）;
flag_2=（（sqrt（（X-N/2）.^2+（Y-N/2）.^2））flag_3=（X>（N/2-15））&（X<（N/2+15））;
flag_4=（Y>（N/2-15））&（Y<（N/2+15））;

flag=flag_1|flag_2;%|flag_3|flag_4;

flag1=not（flag）;

?i_P=i.*flag1;
?
?i1=uint8（255/2*i_P）;
?figureimagesc（i1）;
?colormap?gray
?axis?square

?%%%%%%%%%%%%%%%%%%%%%%%%%?三?求出載頻值??%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
?
?I_P=ft2（i_Pdelta）;
?
?[j?k]=find（abs（I_P）==max（max（abs（I_P））））;?
?s0=j（1）;??t0=k（1）;???????????????????????%?找出零級(jí)頻譜坐標(biāo)
?
?[j?k]=find（abs（I_P（1:（N/2-N/128）1:（N/2-N/128）））...
?==max（max（abs（I_P（1:（N/2-N/128）1:（N/2-N/128））））））;??%?找出負(fù)一級(jí)頻譜坐標(biāo)
?s1=j（1）;t1=k（1）;

fx0=s0-s1;???????????????????????????????????%?求出載頻值??????
fy0=t0-t1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%?四?條紋延拓?%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

????i_E=（i-bg）.*flag1;
%???i1=uint8（255*abs（i_E））;
%???figureimagesc（i1）;
%???colormap?gray
%???axis?square

for?m=1:10
?
????I_E=ft2（i_Edelta）;
%??
%??I1=uint8（255*abs（I_E））;
%??figureimagesc（I1）;
%??colormap?gray
%??axis?square
%??
?%%%%%%%%%%%%%%%%%%%?找出一級(jí)和負(fù)一頻譜坐標(biāo)?%%%%%%%%%%%%%%%%%%%%%%
?
[j?k]=find（abs（I_E（1:（N/2-N/128）1:（N/2-N/128）））...
==max（max（abs（I_E（1:（N/2-N/128）1:（N/2-N/128））））））;??%?在第三象限內(nèi)找出一級(jí)頻譜坐標(biāo)
s11=j（1）;??t11=k（1）;?

[j?k]=find（abs（I_E（（N/2+N/128）:N（N/2+N/128）:N））...
==max（max（abs（I_E（（N/2+N/128）:N（N/2+N/128）:N）））））;??%?在第二象限內(nèi)找出負(fù)一級(jí)頻譜坐標(biāo)
s22=j（1）+N/2+N/128-1;??t22=k（1）+N/2+N/128-1;?

%%%%%%%%%%%%%%%%%%%%%??濾波器設(shè)計(jì)?%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

r1=1/16;
H1=（（sqrt（（X-s11）.^2+（Y-t11）.^2））

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%??濾波??%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I_E=I_E.*H1;

?%%%%%%%%%%%%%%%%%%%%%%%%%?反傅里葉變換?%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
?
i_E=ift2（I_Edelta_f）;


?i_E=i_E.*flag+（i-bg）.*flag1;
?
end
?
?i1=uint8（255*i_E）;
?figureimagesc（i1）;
?colormap?gray
?axis?square

?屬性????????????大小?????日期????時(shí)間???名稱(chēng)
-----------?---------??----------?-----??----

?????文件???????2668??2012-10-17?10:48??explam-fft.m

-----------?---------??----------?-----??----

?????????????????2668????????????????????1