求信号的互相关函数时遇到的问题

Minghanzhang

各们高人:
        我在求两信号的互相关函数时出了问题,程序如下,希望各位高手指点迷津!
         [Pxy1,f] = csd(x1,y1,8192,8000,256,128);
         Rxy1 = ifft(Pxy1);
其中X1和Y1表示两个长度相同的时域信号,我想求他们的互相关函数,我采用间接法
来求,即先对信号X1和Y1用CSD命令来求得互功率谱密度,再经过IFFT来求出它们的
互相关函数.但我求出的互相关函数有复数值,如何求出与XCORR函数求出的类型一
样,也就是说是实数域中的.
我可能对这两个函数的理解不深,肯请高人指教,不慎感激!!

Minghanzhang

另外,我运用了以下方式求互相关函数X(k)=fft(x1(n));
Y(k)=fft(x2(n));
Pxy = conj(X(k)).*Y(k);
Rxy=ifft(Pxy);

不知道对不对,高人们帮我判断一下,谢谢!!!

songzy41

原帖由 Minghanzhang 于 2008-8-20 18:32 发表
各们高人:
        我在求两信号的互相关函数时出了问题,程序如下,希望各位高手指点迷津!
         [Pxy1,f] = csd(x1,y1,8192,8000,256,128);
         Rxy1 = ifft(Pxy1);
其中X1和Y1表示两个长度相同的时域信 ...

用下语句的csd，实际上是求一个平均的互功率谱，每帧长256，2帧之间重叠128，每帧以8192进行FFT变换。经csd得到的Pxy1是在0-4000之间（即正频率部分）的复数谱值。
[Pxy1,f] = csd(x1,y1,8192,8000,256,128);
Pxy1是复数，又没有满足共轭对称的条件，经IFFT变换后，当然是复数：
Rxy1 = ifft(Pxy1);
所以这样的结果与xcorr的结果相差甚远。

songzy41

原帖由 Minghanzhang 于 2008-8-20 18:38 发表
另外,我运用了以下方式求互相关函数X(k)=fft(x1(n));
Y(k)=fft(x2(n));
Pxy = conj(X(k)).*Y(k);
Rxy=ifft(Pxy);

不知道对不对,高人们帮我判断一下,谢谢!!!

这个结果同样不对，这计算的结果是循环互相关（circular correlation），与xcorr的线性互相关不同。要想得到与xcorr一样的结果，可以这样设置：
N=max(length(x1),length(x2));
X(k)=fft(x1,2*N-1);
Y(k)=fft(x2,2*N-1);
Pxy = conj(X(k)).*Y(k);
Rxy=real(ifft(Pxy));

Minghanzhang

感谢songzy41的指点!谢谢,兄弟我不慎感激!

Minghanzhang

高人,我按你的程序运行了下,结果发现两者运行结果不一样,你给的程序运行结果是峰值出现在两头,也就是N=0各N=N的位置;
直接利用XCORR(X,Y)得出的是峰值出现在N/2处.还有一处不解的是你给的程序中Rxy=real(ifft(pxy)),只示实部,那虚部呢?兄弟
我初学,希望高人能指点一二,谢谢!!!

antonylau

原帖由 songzy41 于 2008-8-20 20:36 发表

这个结果同样不对，这计算的结果是循环互相关（circular correlation），与xcorr的线性互相关不同。要想得到与xcorr一样的结果，可以这样设置：
N=max(length(x1),length(x2));
X(k)=fft(x1,2*N-1);
Y(k)=fft(x ...

这个好像是信号分别在前后补零，然后共轭相乘，再ifft
以前有个帖子讲到

songzy41

原帖由 antonylau 于 2008-8-21 10:40 发表


这个好像是信号分别在前后补零，然后共轭相乘，再ifft

antonylau 说得对的，应把x和y分别是前后补零，然后再处理。以下附上程序、数据和计算结果，可看到xcorr计算得的（红色线）和FFT计算得的（蓝色线）完全重合在一起。

x=load('data1.txt');
y=load('data2.txt');
N=length(x);
n=1:N;
%subplot 211; plot(n,x);grid;
%subplot 212; plot(n,y);grid;
R=xcorr(x,y,'biased');
%figure
m=-N+1:N-1;
plot(m,N*R,'r','linewidth',2); hold on;
y=[zeros(1,N-1) y'];
X=fft(x',2*N-1);
Y=fft(y,2*N-1);
Pxy=X.*conj(Y);
Rxy=real(ifft(Pxy));
plot(m,Rxy); grid;
legend('xcorr','fft')

Minghanzhang

恩人,太谢谢你了哈!现在总算搞明白这个函数.此外在请教一下高人,我看别人做的仿真是在不同信噪比下得出的结果,
比如测试一个算法对声音信号的识别能力,不同信噪比下应该得出不同的结果.这种信噪比应该如何设置啊?谢谢!!!

[ 本帖最后由 zhangnan3509 于 2008-8-21 18:35 编辑 ]

songzy41

可以用 degrade 函数，增加噪声达到指定的信噪比。

Minghanzhang

谢谢,但是我help degrade后,查无此函数啊!!!呵呵:lol

songzy41

原帖由 Minghanzhang 于 2008-8-21 20:36 发表
谢谢,但是我help degrade后,查无此函数啊!!!呵呵:lol

我给你附上这程序:
% [Y] = degrade(X,SNR,NOISE)
% adds SNR dB NOISE to the speech signal X and returns the noisy signal in Y
% By default SNR=10dB and NOISE is AWGN.
function [Y] = degrade(X,SNR,NOISE)
if isstr(X)==1, X=read(X); end
if nargin < 3, NOISE=randn(size(X)); end
if nargin < 2, SNR=10; end
if length(NOISE) < length(X),
   disp('Length of speech signal X should be greater than noise NOISE');
   break
end
signal_power = 1/length(X)*sum(X.*X);
noise_variance = signal_power / ( 10^(SNR/10) );
Y=X+sqrt(noise_variance)/std(NOISE)*NOISE(1:length(X));
Y=32000/max(abs(Y))*Y;

Minghanzhang

谢谢songzy41这么长时间对我的支持,兄弟由衷的表示感谢!不过小弟又有问题出来了,
你给的这个程序与我直接调用:y = AWGN(x,SNR,'measured')有什么区别,嘿嘿,高人你
有时间的话在解答吧!:handshake

80423989@qq.com

不知道矩阵的相关函数如何利用fft求解呢？而且相关的定义不同，求法也不同。比如你这里的是r(k)=x(k1)*x(k1-k)的这种定义，所以fft是fft(x).*conj(y).我想请教对于定义为类似r(k)=x(k1)*x(k1+k)这种相关，如何利用ft求呢？x,y都是矩阵，是二维的。