相干函数的基本概念及其案例

一、相干函数的基本概念

设有一个系统，传递函数为H，平稳随机过程输入为x（t），输出为y（t），则有y(t)=h(t)*x(t)

其中：h（t）是传递函数H的脉冲响应，’*‘表示卷积。在频域上有Y(f)=H(f)X(f)或Pyy(f)=H^2(f)*Pxx(f)

式中：Pxx（f）是x（t）的自功率谱；Pyy（f）是y（t）的自功率谱。在已知输入x（t）和输出y（t）时，可以求出传递函数H为

 

在离散条件下，上式转换为

上式中频率f用离散频率索引来表示，f=k△f。
系统的相干函数定义为

 

其中：Pxx（f)是x（t)的自功率谱；Pyy（f)是y（t)的自功率谱；Pxy（f)是x（t)和y（t)的互功率谱。相干函数Yxy（f）可用来描述这两个信号在各频率点处的相关程度。

从相干函数Yxy（f）可以确定输出信号y（t）在多大程度上来自于输入信号x（t）。若Yxy（f）=1，说明输出完全来自于输入，且系统必为线性系统；若Yxy（f）<1，对于线性系统表明在频率点f处，输出谱Pyy（f）有多少成分来自于输入谱Pxx（f），其余部分可能来自于另外的信号源或噪声；若，Yxy（f）=0，则x（t)和y（t)完全不相干。

离散条件下，上式可以转换为

与传递函数一样，上式中频率f用离散频率索引来表示，f=k△f。

二、相干函数的估算
名称：mscohere
功能：对两个数据序列估算相干函数
调用格式：
Cxy=mscohere(x,y)
Cxy = mscohere(x,y,window)
Cxy = mscohere(x,y,window,noverlap)
[Cxy,W] = mscohere(x,y,window,noverlap,nfft)
[Cxy,F] = mscohere(x,y,window,noverlap,nfft,fs)
[.……] = mscohere(x,Y,...,'twosided')
mscohere(.….)
说明：输人参数：x和y是两列数据，表示系统的输人和输出；window是窗函数，默认时用海明窗；noverlap是分段时两段之间的重叠部分；nfft是傅里叶变换时的长度；fs是采样频率；一般x和y为实数序列时求出的相干函数是单边谱，若要求双边谱可用参数'twosided'。输出参数：Cxy是计算出的相干函数；W是归一化的角频率矢量，在单边谱时在[0，pi]之间，双边谱在[0，2*pi]之间；F是频率矢量，只有在输入参数中含有fs才给出F，它是实际频率值。当没有输出参数时，将显示出相干函数的图谱。

案例、有2个FIR滤波器，它们的系数分别为h=ones(1,10)/10和h1=fir1（30,0.2,rectwin（31))。由随机数列通过这两个滤波器产生输出序列为x（n)和y（n)，再用x(n)和y（n)以 welch法求自谱及互谱，并计算出相干函数，然后调用mscohere函数计算出x（n)和y（n)的相干函数，比较这两种方法的计算结果。程序如下：

clear all; clc; close all;
 
randn('state',0);                  % 随机数初始化
h = ones(1,10)/10;                 % 滤波器1系数
h1 = fir1(30,0.2,rectwin(31));     % 滤波器2系数
r = randn(16384,1);                % 产生随机数
x = filter(h,1,r);                 % 产生第1路信号x
y = filter(h1,1,r);                % 产生第2路信号y
 
N=length(x);                       % 数据点长度
[H,wh]=freqz(h,1);                 % 滤波器1的响应函数
[H1,wh1]=freqz(h1,1);              % 滤波器2的响应函数
 
wind=hamming(1024);                % 设置海明窗,窗长1024
noverlap=512;                      % 重叠长度
Nfft=1024;                         % FFT变换长度
PY1=pwelch(x,wind,noverlap,Nfft);  % 求第1路信号自谱
PY2=pwelch(y,wind,noverlap,Nfft);  % 求第2路信号自谱
[CY12,w1]=cpsd(x,y,wind,noverlap,Nfft);   % 求第1路和第2路信号的互谱
Co12=abs(CY12).^2./(PY1.*PY2);     % 按(8-5-5)计算相干函数
[CR,w2]=mscohere(x,y,wind,noverlap,Nfft); % 调用mscohere函数计算相干函数
mcof=max(abs(Co12-CR))             % 求两种方法的差值
% 作图
figure(1)
subplot 211; plot(x,'k'); title('第1路信号x波形');
ylabel('幅值'); xlabel('样点'); axis([0 N -1.2 1.2]);
subplot 212; plot(y,'k'); title('第2路信号y波形');
ylabel('幅值'); xlabel('样点'); axis([0 N -1.2 1.2]);
set(gcf,'color','w'); 
figure(2)
subplot 211; plot(wh/pi,20*log10(abs(H)),'k'); grid;
ylim([-60 10]); title('滤波器1幅值响应曲线');
ylabel('幅值/dB'); xlabel('归一化频率/pi');
subplot 212; plot(wh1/pi,20*log10(abs(H1)),'k'); grid;
ylim([-70 10]); title('滤波器2幅值响应曲线');
ylabel('幅值/dB'); xlabel('归一化频率/pi');
set(gcf,'color','w'); 
figure(3)
plot(w1/pi,Co12,'r','linewidth',2);
hold on; grid;
plot(w2/pi,CR,'k');
legend('调用自谱和互谱','调用mscohere',3)
title('两种方法求得相干函数比较');
ylabel('幅值'); xlabel('归一化频率/pi');
set(gcf,'color','w'); 

运行结果如下：

 

 

运行程序先得到x(n)和y(n)两序列的波形图，上图所示；然后又给出了两个FIR滤波器的幅频响应，上图所示。用两种方法计算了相干函数，并把这两个相干函数显示在一张图上，上图所示，可以看出这两条曲线重合得很好。同时计算了两种方法的差值，得到两种方法的差值为0，说明这两种计算方法是一致的。

参考文献：MATLAB数字信号处理85个实用案例精讲——入门到进阶；宋知用（编著）