声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 2172|回复: 2

[综合] 自编SCOT互相关时延估计matlab程序

[复制链接]
发表于 2012-3-22 17:03 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
clear all;clc;close all;
N=1024;%信号长度
fs=1000;%采样频率
n=0:N-1;
t=n/fs;%时间序列

a1=5;%信号幅度
a2=5;%信号幅度
d=26;%延迟点数

x1=a1*sin(2*pi*10*t)+1.9*sin(2*pi*18*t)+2.8*sin(2*pi*55*t);
x1=x1+randn(size(x1));


x2=a2*sin(2*pi*10*(n+d)/fs)+1.9*sin(2*pi*18*(n+d)/fs)+2.8*sin(2*pi*55*(n+d)/fs);
x2=x2+randn(size(x2));
%x2=x2.*hamming(N);

subplot(311);
plot(t,x1,'r');
axis([-0.1 1.5 -6 6]);
%hold on;
subplot(312);
plot(t,x2,'g');
axis([-0.1 1.5 -6 6]);
%legend('x1信号','x2信号');
%xlabel('时间/s');ylabel('x1(t)x2(t)');
%title('原始信号');grid on;
%hold off;

%SCOT互相关算法
X1=fft(x1,2*N-1);
X2=fft(x2,2*N-1);
Sxy=X1.*conj(X2);
%Cxy=fftshift(ifft(Sxy));
%Sx=abs(X1).^2;
%Sy=abs(X2).^2;
%Pxy=Sxy./((sqrt(Sx.*Sy))+eps);
Pxy=Sxy./(abs(X1).*abs(X2)+eps);
Cxy=fftshift(real(ifft(Pxy));

t1=(0:2*N-2)/fs;
subplot(313);
plot(t1,Cxy,'b');
title('互相关函数');
xlabel('时间/s');ylabel('Cxy(t)');grid on;

[max,location]=max(Cxy)%求出最大值及其所在位置
D=location-N
Delay=D/fs

基于SCOT的互相关时延估计,用这个程序仿真得到的估计值每次都不一样,  请大侠看看是怎么回事。

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

发表于 2012-8-3 19:21 | 显示全部楼层
回复 1 # lifeman001 的帖子

这是我参照别人的代码改的一段代码,看对你是否有帮助
clear;

N=1024;%信号长度
fs=1000;%采样频率
n=0:N-1;
t=n/fs;%时间序列
ts = 1/fs * (-N + (1 : 2*N - 1)); %互相关时间序列

a1=5;%信号幅度
a2=5;%信号幅度
d=26;%延迟点数

x1=a1*sin(2*pi*10*t)+1.9*sin(2*pi*18*t)+2.8*sin(2*pi*55*t);
x2=a2*sin(2*pi*10*(n+d)/fs)+1.9*sin(2*pi*18*(n+d)/fs)+2.8*sin(2*pi*55*(n+d)/fs);

%x2=awgn(x1./4,-3);                %噪声强度大于信号
%x2=x2 .* hamming(N);

x=awgn(x1,20);                     %加入噪声
y=awgn(x2,0.001);

S1=fft(x,2*N-1);
S2=fft(y,2*N-1);
X = S1.* conj(S2);
R1 =real(fftshift(ifft(X)));

subplot(211),plot(ts,R1);
xlabel('时间/s');ylabel('R1(t)');
title('ifft计算结果')

Cc=xcorr(x(1:N),y(1:N));                     
subplot(212),plot(ts, Cc);
xlabel('时间/s');ylabel('Cc(t)');
title('xcorr(x(1:N),y(1:N))');

c1=max(abs(R1));
a1=find(abs(R1)==c1) %找出最大点
m1=(a1-length(R1)/2)*1/fs %换算出信号的延迟时间

c2=max(abs(Cc));
a2=find(abs(Cc)==c2) %找出最大点
m2=(a2-length(Cc)/2)*1/fs %换算出信号的延迟时间
发表于 2012-8-29 16:44 | 显示全部楼层
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-5-2 18:47 , Processed in 0.090690 second(s), 20 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表