|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
本帖最后由 wdhd 于 2016-9-12 14:16 编辑
对于等间距的一系列麦克风,同时采集一组数据。由于状态随着距离在传播,这一系列麦克风采集到的数据具有相关性。只对其相关性分析,可以得到随距离和时间变化时的相关系数等值线分布,相关系数图中可以也看到有两个斜率存在,距离/时间可以得到相应的速度。
现在想利用公式计算得到这个速度,U/X=w/phase;
于是想根据互谱求得相位差,根据距离的关系得到速度。最主要的是求得相位/距离这个斜率。
1.直接对数据进行分析时发现相位跳跃厉害,于是先自己生成信号,看斜率。
2.代码为自己生成一系列信号,求得互谱相位差,画出相位差-距离的图像,发现在不跳跃的地方线性度还可以。应如何改进这个跳跃呢?
dx=0.01;
fs=1000;
x=dx:dx:0.5;
dt=1/1000;
t=0:1/1000:5;
c=20;
%==============================signal generation==========================
for i=1:length(x)
signal(:,i)=sin(2*pi*100*(t-x(i)/c))+rand(1,length(t));
signal1=sin(2*pi*100*t)+rand(1,length(t));
nfftcpsd=2^nextpow2(2000);
ff=0:fs/nfftcpsd:fs/2-fs/nfftcpsd;
%==============direct cpsd method========================================
[Pxy,W]=cpsd(signal1(2000:4000),signal(2000:4000,i),hann(2000),[],nfftcpsd);
Pxy_absdir(:,i)=abs(Pxy);
%find the phase
Pxy_phasedir(:,i)=atan(imag(Pxy)./real(Pxy));
% Pxy_phase(:,num)=angle(Pxy);
Pxy_phaseunwrapdir(:,i)=unwrap(Pxy_phasedir(:,i));
end
%====================crosscorrelation and then fft=======================
%========================correlation channel1 and channelother=============
y1=[];
y2=[];
xcorrelation=[];
nlags=1000;
for j=1:length(x)
y1=signal1(2000:4000);
y2=signal(2000:4000,j);
temp=crosscorr(y1,y2,nlags);
xcorrelation(:,j)=temp;
end
%=========================fft================================================
nfft=2^nextpow2(8000);
f=0:fs/nfft:fs/2-fs/nfft;
%=============find the index for two specific frequency=============
[fmin1,index1]=min(abs(100-f));
[fmin2,index2]=min(abs(170-f));
%===================================================================
for num=1:length(x);
tempxcorr=xcorrelation(:,num);
fftcorr=fft(tempxcorr,nfft);
Pxy_abs(:,num)=abs(fftcorr);
%find the phase
Pxy_phase(:,num)=atan(imag(fftcorr)./real(fftcorr));
% Pxy_phase(:,num)=phase(fftcorr);
Pxy_phaseunwrap(:,num)=unwrap(Pxy_phase(:,num));
Phase1result(num)=-Pxy_phaseunwrap(index1,num);
Phase2resultunwrap(num)=-Pxy_phaseunwrap(index2,num);
Phase2result(num)=Pxy_phase(index2,num);
end
%=====plot 100Hz==============
plot(x,Phase1result)
hold on;
plot(x,Pxy_phaseunwrapdir(206,:),'r')
[ 本帖最后由 antonylau 于 2008-8-21 09:59 编辑 ] |
-
相关分析图
-
相位图比较
|