相关函数的一些应用，欢迎跟帖讨论

花如月

相关函数的应用
  相关函数应用很广，例如噪声中信号的检测、信号时延长度的检测、隐含周期性的检测、信号相关性检验等。而在平稳随机信号里应用就更多了（比如估计信号的功率谱）。由于水平有限，这里举2个简单的应用例子，权作抛砖引玉。希望大家多提意见，积极参与讨论
一 相位延迟的计算
问题描述：
设x(n)和y(n)是2个传感器上得到的观测序列,即:
x(n)=s(n)+w1(n)     (1)
y(n)=s(n-D)+w2(n)   (2)
其中,s(n)是未知信号,s(n-D)为s(n)的移位形式,w1(n)和w2(n)是未知观测噪声。
时延估计问题是通过有限个观测值{x(n),y(n),1≤n≤N}来估计时延参数D。
现有的估计办法：
（1）互相关函数的方法（适用于观测噪声w1(n)与w2(n)为零均值噪声,它们彼此独立,并且与s(n)独立。）
估计的原理：对x(n)和y(n)求互相关，互相关峰值时刻对应的t即为时延D。
（2）基于LMS自适应算法的瞬态信号时延估计
（3）离散小波变换或者小波包的方法
（4）基于人工神经网络的信号时延估计方法
（5）基于小波变换的模极大值方法的时延估计
后边几种方法比较麻烦请查阅相关的文献
各种方法的比较
相关函数法时延估计,方法简单,易于实现,但在背景噪声超过实际的声波发生信号强度时,会出现比较大的测量误差;基于人工神经网络的时延估计,根据样本数据训练神经网络,训练完毕后,可实现含噪信号到纯有用信号的非线性映射,然而一旦背景噪声的类型发生改变时,神经网络需
要重新进行训练,同时神经网络训练方法的选择!神经网络结构的最优化等问题还有待于进一步的探讨;基于小波变换的信号时延估计,目前已经应用在各个领域中,从其实际的应用效果来看,是一种值得推广使用的信号时延估计方法,相比于其他两种方法,小波变换技术更适合于强噪声环境下的声波信号的时延估计"（参考文献来自中国知网）

在信号没有噪声的情况下，时延的估计其实就是相位的估计。这个问题版友songzy41 在http://forum.vibunion.com/thread-39455-1-1.html里了解答。下边是他的程序
n=99;                              %设定每周期数据采集点数
T=6;                               %采样周期数
t=0:2*pi/(n-1):2*T*pi;             %采样数
N=length(t);
fs=1/(n-1);
y1=4*sin(t);                       %信号1
y2=8*sin(t+pi/6);                  %信号2，相位差取pi/6
Cc=xcorr(y1,y2);                   %求互相关函数
plot(Cc);
[y_max,t_max]=max(Cc)              %找出Cc的最大值及对应的t_max
sample_delay=t_max-N;              %计算与中心点相差的样点数
T_sig=(N-1)/6;                     %求信号一个周期的样点数,一个周期相对应于2*pi
delay1=pi/6                        %pi/6的弧度
delay2=2*pi*sample_delay/T_sig     %计算与中心点相差的样点数所对应的弧度值

看到一个版友问到用FFT 的时延估计问题
“我看到了一个关于“用FFT快速实现时延估计”的程序，如下：设两路信号分别为x,y且
S1=fft(x,2048) , S2=fft(y,2048)
X = S1.* conj(S2);
Y = abs(S1) .* abs(S2);
Y = Y+ eps;
G = X ./ Y;
R = ifft(G);
现在求出R后，还需怎样计算才能得出两信号的时延呀？”
这个程序估计时延的理论基础其实还是相关函数估计时延。而相关函数的计算利用的的是时域相关定理。
下面验证二种方法是一致的，也就是验证时域相关定理。
clear
n=99;                              %设定每周期数据采集点数
T=6;                               %采样周期数
t=0:2*pi/(n-1):2*T*pi;             %采样数
N=length(t);
fs=1/(n-1);
x1=4*sin(t);                       %信号1
y1=4*sin(t+pi/6);                  %信号2，相位差取pi/6
x=awgn(x1,20);                     %加入噪声
y=awgn(y1,0.001);
x2=awgn(x1./4,-3);                %噪声强度大于信号
subplot(311),plot(t,x)
subplot(312),plot(t,y)
tic
Cc=xcorr(x,y);                      %直接求互相关函数
toc
Cc1=xcorr(x1,y1);  
subplot(313),plot(Cc,'r.'),hold on
plot(Cc1,'-');
tic                                 %根据时域相关定理求互相关
S1=fft(x,2*N-1);
S2=fft(y,2*N-1);
X = S1.* conj(S2);
R1 =real(fftshift(ifft(X)));
toc
figure(2)
subplot(211),plot(Cc),title('xcorr结果')
subplot(212),plot(R1);
title('ifft计算结果')
figure(3)
plot(xcorr(x2,'unbiased'));            
由figure2可以看到2者计算的结果是吻合的，程序调试过程中得到songzy41老师的指导。再次表示感谢。
我们知道wgn是调用randn产生白噪声的，而randn产生的是均值为零，方差为1的高斯白噪声。所以上边的程序检测时延是比较准确的，即使在信号被噪声淹没的情况下，相关函数的仍然可以检测到。可以通过将加入噪声的信噪比设置为负值，然后观察效果。通过统计的时间来看，直接用xcorr求互相关比用时域相关定理法求快一些。
Elapsed time is 0.001933 seconds.
Elapsed time is 0.004124 seconds.

二 噪声中信号的检测
原理：对一噪声污染严重的信号，做自相关。从figure3可以辨别出原始信号包含正弦信号。这是一个《数字信号处理》 胡广书编 清华大学出版社 38页的例子，不过图和书上给的有些出入。书上在0点取得最大值，后边就比较正常了，而这个画出来的图像是对称的，并且在600左右一个尖峰。不知道是怎么回事，由于小弟信号处理的基础非常差，所以犯错是再所难免的，望各位老师、前辈版友不惜赐教。

欢迎大家批评指正，踊跃跟帖！

相关函数概念和计算请看eight8版主的这个帖子
[原创]自相关函数和互相关函数计算和作图的整理及一点心得

[ 本帖最后由 花如月 于 2007-7-27 17:02 编辑 ]

w89986581

“所以上边的程序检测时延是比较准确的，即使在信号被噪声淹没的情况下，相关函数的仍然可以检测到。”是因为构造的干扰是高斯白噪声。有些情况下，噪声是具备时间、空间相关性的，检测起来就比较麻烦了，呵呵。

weipeng

首先、非常感谢您的细致讲解。
其次、“y1=4*sin(t+pi/6);    %信号2，相位差取pi/6”
中pi/6是任取的还是默认的？
第三、“subplot(313),plot(Cc,'r.'),hold on
plot(Cc1,'-');”中Cc与Cc1的图形完全重合，这说明什么呢？是对加入噪声与不加入噪声的信号，做互相关的结果是相同吗？

花如月

第一个问题：相位差是自己定的，最终程序会求出这个相位差。
第二个问题：你的理解是对的

weipeng

这两个图为什么不一样呢？
[x,fs,NBITS]=wavread('d:\ttt.wav');
[y,fs,NBITS]=wavread('d:\ttt-1.5.wav');
S1=fft(x,2048);
S2=fft(y,2048);
X = S1.* conj(S2);
R1 =real(fftshift(ifft(X)));
Cc=xcorr(x,y);                     
subplot(211),plot(R1);
title('ifft计算结果')
subplot(212),plot(Cc);
title('xcorr结果');

花如月

对照数据自己检查一下，注意保证2二者的长度相等。应该是个小问题，自己排查一下

weipeng

所谓的长度相等是指两个声音文件的时间长度相等还是N=length(　)的长度呀？我用手动录制的声音很难长度相等，是不是还得用软件处理一下呀？
我又实验了一下，把同一个声音复制了一份ttt.wav和ttt_copy.wav,程序如下：
[x,fs,NBITS]=wavread('d:\毕业设计\实验\ttt.wav');
[y,fs,NBITS]=wavread('d:\毕业设计\实验\ttt_copy.wav');

S1=fft(x,2048);
S2=fft(y,2048);
X = S1.* conj(S2);
R1 =real(fftshift(ifft(X)));
Cc=xcorr(x,y);                     
subplot(211),plot(R1);
title('ifft计算结果')
subplot(212),plot(Cc);
title('xcorr结果');
这两个图还是不太一样，哪错了呢？请指教。

花如月

第一个数据序列只有2000点左右，第二个都60多万点了。当然会不一样了，可以考虑数据去短些。让fft的点数等于互相关运算后的序列的长度，才具有可比性

[ 本帖最后由 花如月 于 2007-7-26 21:58 编辑 ]

gladstone

那里有自相关函数消除白噪音的证明啊？
含有白噪音的锯齿序列经自相关运算后呈现正弦样式是怎么回事啊？
如能赐教，不胜感谢。

花如月

（1）自相关函数是不能消除噪声的，他只是可以检测隐含在噪声中的信号以及周期性
（2）具体原因我也说不上来，不过觉得是正常的。你可以把你的程序发上来，让大家帮你分析一下。

weipeng

您好，在您和一位朋友的指点下，我的程序有点进展了。现在两个图的大致形状相同，但Y轴的幅值一样。程序如下：
[x,fs,NBITS]=wavread('d:\毕业设计\实验\ttt.wav');
[y,fs,NBITS]=wavread('d:\毕业设计\实验\ttt_copy.wav');
S1=fft(x,2048);
S2=fft(y,2048);
X = S1.* conj(S2);
R1 =real(fftshift(ifft(X)));
Cc=xcorr(x(1:2000),y(1:2000));                     
subplot(211),plot(R1);
title('ifft计算结果')
subplot(212),plot(Cc);
title('xcorr结果');

还需怎么改进呢？谢谢指教！

花如月

（1）如果想要完全形同必须是对同一个序列
（2）在绘制结果的时候保证序列长度一致
幅度不一样确实有些奇怪，程序应该没有问题，那个例子是完全吻合的。你再仔细检查一下
我电脑这会儿matlab跑不动了。空了我找个wav试试

[ 本帖最后由 花如月 于 2007-7-30 16:57 编辑 ]

weipeng

您　好，
Cc=xcorr(x(1:2048),y(1:2048));   % －－－－1         
dd=xcorr(x,y,2048);%－－－－－－2

这两句是不是作用和效果一样呀？
用2画出的图形状与上面的图差不多，但Y轴辐值是－500到500。这是怎么回事呢？

对了，求两信号的时间差（时延），您说的是“对x(n)和y(n)求互相关，互相关峰值时刻对应的t即为时延D。”一位版友这样写的
c1=max(abs(R1));
a1=find(abs(R1)==c1) %找出最大点
m1=(a1-length(R1)/2)*1/fs %换算出信号的延迟时间
这样写对吗？我感觉好象不太对，请指教。

花如月

你的问题看一下help完全是可以解决的，最后换算信号的延迟时间也是对的。我的matlab依然有问题，所以还是没法帮你试

weipeng

您好，经常问这种简单问题，实在不好意思，但为了学点知识，还得问：
我对一个声音信号和它的复制文件（两个一样的信号）进行计算时延实验，两个信号的时延是不是应该是0呀？
但我用如下的程序，三种方法进行比较，得出来的都不是0，不知为什么？
clear;
[x,fs,NBITS]=wavread('d:\毕业设计\实验\ttt.wav');
[y,fs,NBITS]=wavread('d:\毕业设计\实验\ttt_copy.wav');
S1=fft(x,4096);
S2=fft(y,4096);
X = S1.* conj(S2);
R1 =real(fftshift(ifft(X)));
Cc=xcorr(x(1:2048),y(1:2048));                     
dd=xcorr(x,y,2048);
subplot(311),plot(R1);
title('ifft计算结果')
subplot(312),plot(Cc);
title('xcorr结果');
subplot(313),plot(dd);
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 %换算出信号的延迟时间

c3=max(abs(dd));
a3=find(abs(dd)==c3) %找出最大点
m3=(a3-length(dd)/2)*1/fs %换算出信号的延迟时间
三种方法找出的最大点分别是a1=2049; a2=2048; a3=2049;
得出的时延分别是：m1= 4.5351e-005; m2=2.2676e-005; m3=2.2676e-005
这是怎么回事呢？