|
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % two_dimension_holospectra_2006_RIDC_WFJ
- % 输入:data_x为H方向信号;data_y为V方向信号;
- % 输入:Fs为采样频率;an为H传感器与X轴夹角;
- % 输入时请注意H方向信号为两传感器中与X轴夹角较小的信号;
- % V方向信号为两传感器中与X轴夹角较大的信号;
- % 若按观察传感器的方向观察轴心轨迹请将direct置0;若要从反方向观察请将direct置1;
- % inter_freq为插值后的工频频率.
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- function [inter_freq]=holo2D('PIP2.DH','PIP2.DV',2000,45,0)
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- H=textread (data_x,'%f',1024);V=textread (data_y,'%f',1024);H=H';V=V';
- if length(H)<1024
- H(length(H)+1:1024)=0;
- end
- if length(V)<1024
- V(length(V)+1:1024)=0;
- end
- mean_H=mean(H);mean_V=mean(V);Hmean=H-mean_H;Vmean=V-mean_V;rHmean=Hmean;rVmean=Vmean;
- i1=1:512;w1=0.5*(1-cos(2*pi*i1/(1024+1)));i2=513:1024;w2=w1(1024-i2+1);w=[w1,w2];
- Hmean=Hmean.*w;Vmean=Vmean.*w;Hmean=Hmean'-mean(Hmean);Vmean=Vmean'-mean(Vmean);
- t=(0:1023)./Fs;
- figure(1)
- subplot(2,1,1);plot(t,rHmean);
- grid on;set(gcf,'color',[1 1 1]);xlabel('时间 t');ylabel('H 幅值 \mum');title('时域信号');
- axis([0 1023/Fs min(rHmean) max(rHmean)]);
- % set(gca,'gridLineStyle','-','FontWeight','Bold')
- subplot(2,1,2);plot(t,rVmean);
- grid on;set(gcf,'color',[1 1 1]);xlabel('时间 t');ylabel('V 幅值 \mum');
- axis([0 1023/Fs min(rVmean) max(rVmean)]);
- % set(gca,'gridLineStyle','-','FontWeight','Bold')
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- const=pi/180;deltf=Fs/1024;fs=(0:1023).*deltf;
- yH=fft(Hmean,1024);amplH=2*abs(yH)./512;phH=angle(yH)./const;
- yV=fft(Vmean,1024);amplV=2*abs(yV)./512;phV=angle(yV)./const;
- figure(2)
- subplot(2,1,1);plot(fs(1:512),amplH(1:512));
- grid off;set(gcf,'color',[1 1 1]);xlabel('频率 f');ylabel('H 幅值 \mum');title('频域信号');
- axis([0 fs(512) 0 max(amplH)]);
- % set(gca,'gridLineStyle','-','FontWeight','Bold')
- subplot(2,1,2);plot(fs(1:512),amplV(1:512));
- grid off;set(gcf,'color',[1 1 1]);xlabel('频率 f');ylabel('V 幅值 \mum');
- axis([0 fs(512) 0 max(amplV)]);
- % set(gca,'gridLineStyle','-','FontWeight','Bold')
- [frespeed,amplampl]=ginput(1);
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- fre=frespeed;num_fre=round(fre/deltf)+1;
- [fremax,fren]=max(amplH((num_fre-5):(num_fre+5)));
- fren=fren+(num_fre-5)-1;frequence=fs(fren);
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- sum_ampl=amplH.^2+amplV.^2;
- for i=1:4
- freq=i.*frequence;num_freq=round(freq/deltf)+1;
- [hvmax,hvn]=max(sum_ampl((num_freq-2):(num_freq+2)));hvn=hvn+(num_freq-2)-1;
- H1=amplH(hvn-1);H3=amplH(hvn+1);H2=amplH(hvn);Hfl=fs(hvn-1);Hfr=fs(hvn+1);Hfm=fs(hvn);HI=hvn;angle_yH=angle(yH(HI));
- V1=amplV(hvn-1);V3=amplV(hvn+1);V2=amplV(hvn);Vfl=fs(hvn-1);Vfr=fs(hvn+1);Vfm=fs(hvn);VI=hvn;angle_yV=angle(yV(VI));
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- [Hfa,HAa,HPa,Vfa,VAa,VPa]=hanning_window(H1,H3,H2,Hfm,V1,V3,V2,Vfm,deltf,angle_yH,angle_yV);
- HFa(i)=Hfa;HAmpl(i)=HAa;HPhase(i)=-HPa./const;
- VFa(i)=Vfa;VAmpl(i)=VAa;VPhase(i)=-VPa./const;
- end
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- if (round(frequence/deltf)-5)<5
- frequence=10.*deltf;
- end
- for i=2:round(frequence/deltf)-5
- if sum_ampl(i-1)<SUM_AMPL(I) sum_ampl(i+1)<sum_ampl(i)
- record(i)=sum_ampl(i);
- else
- record(i)=0.001;
- end
- end
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- [max1,no1]=max(record);record(no1)=0;[max2,no2]=max(record);record(no2)=0;
- [max3,no3]=max(record);record(no3)=0;[max4,no4]=max(record);record(no4)=0;
- maxmax=[max1,max2,max3,max4];nono=[no1,no2,no3,no4];
- [min11,min1]=min(nono);minm(1)=maxmax(min1);non(1)=nono(min1);nono(min1)=round(frequence/deltf)-5;
- [min22,min2]=min(nono);minm(2)=maxmax(min2);non(2)=nono(min2);nono(min2)=round(frequence/deltf)-5;
- [min33,min3]=min(nono);minm(3)=maxmax(min3);non(3)=nono(min3);nono(min3)=round(frequence/deltf)-5;
- [min44,min4]=min(nono);minm(4)=maxmax(min4);non(4)=nono(min4);nono(min4)=round(frequence/deltf)-5;
- for i=1:4
- H1=amplH(non(i)-1);H3=amplH(non(i)+1);H2=amplH(non(i));Hfl=fs(non(i)-1);Hfr=fs(non(i)+1);Hfm=fs(non(i));HI=non(i);angle_yH=angle(yH(HI));
- V1=amplV(non(i)-1);V3=amplV(non(i)+1);V2=amplV(non(i));Vfl=fs(non(i)-1);Vfr=fs(non(i)+1);Vfm=fs(non(i));VI=non(i);angle_yV=angle(yV(VI));
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- [Hfa,HAa,HPa,Vfa,VAa,VPa]=hanning_window(H1,H3,H2,Hfm,V1,V3,V2,Vfm,deltf,angle_yH,angle_yV);
- HFa(i+4)=Hfa;HAmpl(i+4)=HAa;HPhase(i+4)=-HPa./const;
- VFa(i+4)=Vfa;VAmpl(i+4)=VAa;VPhase(i+4)=-VPa./const;
- end
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- inter_freq=HFa(1);
- msgmsg=['您选取的工频为',num2str(inter_freq),'Hz,关闭对话框后按任意键继续程序运行。'];
- msgbox(msgmsg,'Prompt','help');
- pause;
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- Ampl_A=HAmpl/2;Ampl_B=VAmpl/2;Angle_A=90-HPhase;Angle_B=90-VPhase;
- Ampl=[Ampl_A,Ampl_B];max_ampl=1.1*max(Ampl);max_amp=fix(max_ampl/10)*10;
- if max_amp==0
- max_amp=8;
- end
- for i=1:length(Ampl_A)
- t=1:361;
- X(i,:)=Ampl_A(i).*sin((t-1).*const+(Angle_A(i))*const);
- Y(i,:)=Ampl_B(i).*sin((t-1).*const+(Angle_B(i))*const);
- x(i,:)=X(i,:).*cos(an.*const)-Y(i,:).*sin(an.*const);
- y(i,:)=X(i,:).*sin(an.*const)+Y(i,:).*cos(an.*const);
- tan_xy(i,:)=y(i,:)./x(i,:);
- end
- for i=1:length(Ampl_A)
- p(i)=Ampl_A(i).^2+Ampl_B(i).^2;q(i)=2*Ampl_A(i)*Ampl_B(i)*sin((Angle_A(i)-Angle_B(i)).*const);
- long_a(i)=(sqrt(p(i)+q(i))+sqrt(p(i)-q(i)))/2;short_b(i)=(sqrt(p(i)+q(i))-sqrt(p(i)-q(i)))/2;
- e(i)=sqrt(long_a(i).^2-short_b(i).^2)./long_a(i);
- Sx(i)=Ampl_A(i)*cos(Angle_A(i).*const)*cos(an.*const)-Ampl_B(i)*cos(Angle_B(i).*const)*sin(an.*const);
- Sy(i)=Ampl_B(i)*cos(Angle_B(i).*const)*cos(an.*const)+Ampl_A(i)*cos(Angle_A(i).*const)*sin(an.*const);
- Cx(i)=Ampl_A(i)*sin(Angle_A(i).*const)*cos(an.*const)-Ampl_B(i)*sin(Angle_B(i).*const)*sin(an.*const);
- Cy(i)=Ampl_B(i)*sin(Angle_B(i).*const)*cos(an.*const)+Ampl_A(i)*sin(Angle_A(i).*const)*sin(an.*const);
- if direct==1
- Sx1(i)=-Ampl_A(i)*cos(Angle_A(i).*const)*cos(an.*const)+Ampl_B(i)*cos(Angle_B(i).*const)*sin(an.*const);
- Cx1(i)=-Ampl_A(i)*sin(Angle_A(i).*const)*cos(an.*const)+Ampl_B(i)*sin(Angle_B(i).*const)*sin(an.*const);
- Sy1(i)=Sy(i);Cy1(i)=Cy(i);
- else
- Sx1(i)=Sx(i);Sy1(i)=Sy(i);Cx1(i)=Cx(i);Cy1(i)=Cy(i);
- end
- SC=(2*Cx(i)*Cy(i)+2*Sx(i)*Sy(i))./(Cx(i).^2+Sx(i).^2-Cy(i).^2-Sy(i).^2);
- angle_d(i)=0.5*atan(SC)./const;tan_angle(i)=tan(angle_d(i).*const);
- cmp=abs(tan_xy(i,:)-tan_angle(i));[mincmp,mino]=min(cmp);
- axis_cmp=x(i,mino).^2+y(i,mino).^2;acmp=abs(long_a(i).^2-axis_cmp);bcmp=abs(short_b(i).^2-axis_cmp);
- if acmp>bcmp
- angle_d(i)=angle_d(i)+90;
- end
- if angle_d(i)<0
- angle_d(i)=angle_d(i)+180;
- end
- if angle_d(i)>=180
- angle_d(i)=angle_d(i)-180;
- end
- if direct==1
- angle_d(i)=180-angle_d(i);
- end
- end
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- if direct==1
- x=-1.*x;
- end
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- cha=[-3 -1 1 3 -3 -1 1 3];
- figure(3)
- subplot(2,1,1);
- hold on
- for i=1:4
- plot(x(i,1)+cha(i).*max_amp,y(i,1),'r.','MarkerSize',10);plot(x(i,:)+cha(i).*max_amp,y(i,:));
- plot([0+cha(i).*max_amp,x(i,1)+cha(i).*max_amp],[0,y(i,1)],'r-');
- %plot([cha(i).*max_amp,cha(i).*max_amp],[-max_ampl,+max_ampl],'k-');
- direc_X1=x(i,1).*cos(-angle_d(i).*const)-y(i,1).*sin(-angle_d(i).*const);
- direc_Y1=x(i,1).*sin(-angle_d(i).*const)+y(i,1).*cos(-angle_d(i).*const);
- direc_X2=x(i,2).*cos(-angle_d(i).*const)-y(i,2).*sin(-angle_d(i).*const);
- direc_Y2=x(i,2).*sin(-angle_d(i).*const)+y(i,2).*cos(-angle_d(i).*const);
- [direc]=dire(direc_X1,direc_Y1,direc_X2,direc_Y2);
- txt=[direc,num2str(i),'X '];
- text(cha(i).*max_amp-0.25*max_amp,-1.2*max_ampl,txt,'FontWeight','Bold');
- end
- axis off;set(gcf,'color',[1 1 1]);axis equal;axis([-4*max_ampl 4*max_ampl -max_ampl max_ampl]);
- plot([-4*max_ampl,4*max_ampl],[-max_ampl,-max_ampl],'k-');plot([-4*max_ampl,4*max_ampl],[max_ampl,max_ampl],'k-');
- plot([-4*max_ampl,-4*max_ampl],[-max_ampl,+max_ampl],'k-');plot([4*max_ampl,4*max_ampl],[-max_ampl,max_ampl],'k-');
- plot([-4*max_ampl,4*max_ampl],[0,0],'k-');
- title('The 2 - Dimensional Holospectrum','FontWeight','Bold');
- hold off
- subplot(2,1,2);
- hold on
- for i=5:8
- plot(x(i,1)+cha(i).*max_amp,y(i,1),'r.','MarkerSize',10);plot(x(i,:)+cha(i).*max_amp,y(i,:));
- plot([0+cha(i).*max_amp,x(i,1)+cha(i).*max_amp],[0,y(i,1)],'r-');
- %plot([cha(i).*max_amp,cha(i).*max_amp],[-max_ampl,+max_ampl],'k-');
- direc_X1=x(i,1).*cos(-angle_d(i).*const)-y(i,1).*sin(-angle_d(i).*const);
- direc_Y1=x(i,1).*sin(-angle_d(i).*const)+y(i,1).*cos(-angle_d(i).*const);
- direc_X2=x(i,2).*cos(-angle_d(i).*const)-y(i,2).*sin(-angle_d(i).*const);
- direc_Y2=x(i,2).*sin(-angle_d(i).*const)+y(i,2).*cos(-angle_d(i).*const);
- [direc]=dire(direc_X1,direc_Y1,direc_X2,direc_Y2);
- txt=[direc,num2str(round(HFa(i)/HFa(1)*100)/100),'X '];
- text(cha(i).*max_amp-0.25*max_amp,-1.2*max_ampl,txt,'FontWeight','Bold');
- end
- axis off;set(gcf,'color',[1 1 1]);axis equal;axis([-4*max_ampl 4*max_ampl -max_ampl max_ampl]);
- plot([-4*max_ampl,4*max_ampl],[-max_ampl,-max_ampl],'k-');plot([-4*max_ampl,4*max_ampl],[max_ampl,max_ampl],'k-');
- plot([-4*max_ampl,-4*max_ampl],[-max_ampl,+max_ampl],'k-');plot([4*max_ampl,4*max_ampl],[-max_ampl,max_ampl],'k-');
- plot([-4*max_ampl,4*max_ampl],[0,0],'k-');
- hold off
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- figure(4)
- hold on
- text(0,30,'序号');text(8,30,'Sx');text(16,30,'Cx');text(24,30,'Sy');text(32,30,'Cy');text(40,30,'离心率');text(48,30,'长轴倾角');
- for i=1:4
- sq=[num2str(5-i),'X'];text(0,3*(i+4),sq);
- Sx_out=num2str(round(Sx1(5-i).*100)./100);text(8,3*(i+4),Sx_out);
- Cx_out=num2str(round(Cx1(5-i).*100)./100);text(16,3*(i+4),Cx_out);
- Sy_out=num2str(round(Sy1(5-i).*100)./100);text(24,3*(i+4),Sy_out);
- Cy_out=num2str(round(Cy1(5-i).*100)./100);text(32,3*(i+4),Cy_out);
- e_out=num2str(round(e(5-i).*1000)./1000);text(40,3*(i+4),e_out);
- angle_out=num2str(round(angle_d(5-i).*10)./10);text(48,3*(i+4),angle_out);
- end
- for i=5:8
- sq=[num2str(round(HFa(13-i)/HFa(1)*100)/100),'X'];text(0,3*(i-5),sq);
- Sx_out=num2str(round(Sx1(13-i).*100)./100);text(8,3*(i-5),Sx_out);
- Cx_out=num2str(round(Cx1(13-i).*100)./100);text(16,3*(i-5),Cx_out);
- Sy_out=num2str(round(Sy1(13-i).*100)./100);text(24,3*(i-5),Sy_out);
- Cy_out=num2str(round(Cy1(13-i).*100)./100);text(32,3*(i-5),Cy_out);
- e_out=num2str(round(e(13-i).*1000)./1000);text(40,3*(i-5),e_out);
- angle_out=num2str(round(angle_d(13-i).*10)./10);text(48,3*(i-5),angle_out);
- end
- axis off;axis equal;axis([0 60 0 33]);box on;set(gcf,'color',[1 1 1]);
- hold off
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- return
复制代码
程序可以这样贴上来,方便别人查看
另外,function [inter_freq]=holo2D('PIP2.DH','PIP2.DV',2000,45,0)
里面参数根本没有传进去
[ 本帖最后由 sigma665 于 2008-5-19 11:03 编辑 ] |
评分
-
1
查看全部评分
-
|