一、相干函数的基本概念
设有一个系统,传递函数为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个实用案例精讲——入门到进阶;宋知用(编著)