【循环谱】基于最小方差无失真响应滤波器组的谱相关密度估计【附MATLAB代码】

文章来源：微信公众号 EW Frontier

摘要

循环平稳信号的谱相关密度（SCD）估计中的先前工作基于已知具有低分辨率和高频泄漏的加权重叠段平均（WOSA）估计器。在本文中，我们解决了估计的SCD使用高分辨率的最小方差无失真响应（MVDR）估计。使用匹配滤波器组（MAFI）的方法来进行谱估计，我们得到SCD估计的基础上经典的WOSA和MVDR。MAFI方法提供了一个共同的基础，允许SCD估计之间的直接比较，并提供了更好的洞察它们之间的关系。数值算例验证了推导的估计量，并定量评估其性能的周期频率分辨率，方差和偏差的SCD估计。引言在一些工程和科学领域中，大多数信号是通过周期性过程（即调制）生成的，因此所生成的信号表现出在时间上具有周期性的（二阶）矩。这些信号被更好地建模为循环平稳的，并且可以通过循环描述符来表征，如循环自相关和谱相关密度[2]。这些循环描述符是基本工具，从中可以获得更复杂的循环描述符，如循环幅度和循环相位，循环相干谱和循环Wigner-Ville谱，这些循环描述符在通信系统，机械振动分析，气象学，经济学，水文学和声学中具有丰富的实际应用[1，5，6，7，8，20]。对于大多数应用来说，在频域中分析循环平稳信号的结构往往更方便和自然。因此，从观测信号的样本中估计谱相关密度（SCD）具有实际意义。文献中提出的大多数SCD估计量都是基于WOSA [9]，众所周知，WOSA [9]具有较差的分辨率特性和较大的频谱泄漏[1] [5]。当我们处理实际系统中常见的信号的短长度观测时，这些限制尤其麻烦。最小方差无失真响应（MVDR）[3，11]在过去十年中在阵列处理和谱估计中受到关注，因为即使在短长度观测下也具有更高的分辨率能力。在这项工作中，我们认为，SCD估计也可以受益于这些高分辨率的属性和评估MVDR作为一个很好的替代SCD估计。

图1理论和估计SCD。

出于评估目的，我们使用匹配滤波器组（MAFI）方法推导出三个估计量，两个基于MVDR，一个基于WOSA [21，16]。MAFI方法提供了一个共同的基础，其中每种方法仅在滤波器设计及其带宽估计方面有所不同。这允许方法之间的直接比较，并简化了它们的实现。在下一节中，我们提出了一个简短的审查SCD，并提出了一个简单的例子，循环平稳信号，稍后用于验证推导的SCD估计。我们还制定了SCD估计问题作为一个联合频谱估计问题（即交叉频谱）使用循环解调。在第4节的MAFI方法之后，我们得到交叉谱估计MVDR和WOSA的基础上，并提出如何使用这些估计估计SCD。在第5节中，我们数值评估这些估计的统计性质的周期频率分辨率，偏差和方差，最后在6中，我们提出了一些讨论的数值结果和未来的研究工作。

谱相关密度

一个零均值随机过程x（n）是循环平稳的（在广义上），如果由下式给出的自相关函数

在时间上是周期性的，具有整数周期N：

由于自相关Rx（n,t）是周期性的，因此它接受傅立叶级数展开，并且可以写为：

其中傅立叶系数：

是循环自相关，α是循环频率。注意，α=0对于循环相关R，其简化为传统的自相关。与传统的谱分析一样，可以表明这些滞后相关系数的傅立叶变换给予谱相关密度（SCD）的提高[5]：

上述公式显示了信号相对于频谱频率和循环频率w的功率分布。在这方面，谱相关密度包含与信号的非平稳特征相关的附加维度。

3.1光谱相关特性

从经典的谱估计理论可知，具有周期N的周期信号的谱也是具有周期的周期序列。因此，光谱相关性的支持区域在w平面中包含在-π≤w=2πk/N≤π之间,-N/2≤k≤N/2。此外，频谱相关密度（5）在双频平面α-w中呈现以下特性[1]：

在双频平面中的这些对称性加上平面中的周期性，将频谱相关的频率支持限制在双频平面的主象限。

3.2循环平稳信号的例子

循环平稳信号的一个简单例子是加性噪声中的谐波正弦波：

其中，假设v(n)具有零均值和零方差的真实的平稳白色噪声是恒定幅度，并且wc和φ分别是（0，π）和（-π，π]中的确定性常数。使用方程（4）和方程（5），我们可以推导出该模型的SCD的表达式为[6，8]：

图2滤波器组对谱密度估计的解释。

3.3循环解调

为了理解这种解释，我们只需将等式（1）给出的自相关函数定义插入等式（4）：

然后定义u(n)和v(n)并且使得：

将它们替换到等式（10）中以获得：

上面的等式示出Rαx(t)是u(n)和v(n)的互相关性，因此从等式（5）得出Sα(w)是u(n)和v(n)的互谱密度。这种解释表明，Sα(w)可以使用任何交叉PSD估计器（如WOSA和MVDR）进行估计。

滤波器组PSD估计

图2中描绘了经由滤波器组的谱估计的概念。基本上我们有一个零均值信号向量

的谱密度Sα(w)并通过一组窄带滤波器

每一个都以不同的频率w操控滤波器的输出由下式给出：

输出功率P（w）为：

其中Rxx是xn的自相关矩阵的估计值。通过测量一个滤波器的输出功率P（w），我们可以通过关系式[12，14]估计频谱密度S（w）：

在该频率处，控制滤波器，其中Bn是有效滤波器带宽的估计。

基于等式（16），然后我们可以通过选择不同的滤波器S（w）和带宽参数来设计不同的估计器。

4.1 WOSA方法

我们可以设计的最简单的滤波器hwc是幅度为1的N样本的矩形窗口：

在一个频率上操纵的该滤波器w变为：

并且具有恒定带宽：

将这些代入等式（16），我们得到：

其中aw称为导向矢量，定义为：

等式（20）是一般的加权重叠频谱平均（WOSA）频谱估计器。在该基本形式中，方程（20）等价于公知的周期图：

通过将信号向量xn分割成长度Q段并对每个段的谱密度进行平均，我们可以获得Bartlet（平均周期图）方法：

其中，M=N/Q是段数，

添加图片注释，不超过 140 字（可选）

我们还可以应用窗口Rxx（滤波器）和/或重叠片段;这些操作将导致其他经典的非参数谱估计方法。在我们的数值示例中，我们将不使用任何预窗口，因为它们会进一步降低WOSA方法已经很差的分辨率，并且将始终使用给出WOSA方法的最小周期泄漏的分段[1]。

4.2 MVDR方法

WOSA使用的平方滤波器具有sinc形状的脉冲响应H（wc），如图3所示，我们可以看到这种形状对P(wc)的估计有直接影响。

图3通过滤波的谱密度估计。

为了提高估计值，MVDR实现了一个窄带滤波器，其脉冲响应在该频率wc处等于1（无失真），同时尽可能降低其他频带处的功率（最小方差）[11]。

该滤波器设计是一个众所周知的优化问题，我们通过最小化功率来解决hwc：

受到hwcHa=1约束，其中a是在（21）中定义的导向矢量。该优化问题具有由[3，11]给出的已知解：

在（15）中替换它给出PSC（功率谱Capon）估计器：

为了得到S（wc），我们必须按照关系式（16）通过滤波器带宽Bn进行P（wc）归一化。滤波器长度的倒数给出了一个简单的Bn估计值，Bn=1/Q.也就是说，这导致了MVDR估计器[3]：

另一个被称为归一化MVDR或NMVDR的估计量是通过使用Bn=aHa导致[13]的关系获得的：

从等式（30）和（31）可以看出，MVDR和NMVDR都是自适应滤波器，因为它们取决于信号特征Rxx。这与WOSA使用的矩形滤波器形成对比，后者是通用的，并且在其构造中不考虑要滤波的信号。4.3互功率谱估计和SCD估计。

图4通过循环迭代法估计SCD。

从滤波的观点来看，根据两个数据向量的互功率谱的估计可以基于两个窄带通滤波器的x1和x2设计并且在相同的频率处操纵。因此，互功率谱可以被推断为滤波器输出在滞后零处的互相关[14]：

其中y1和y2分别为hx1Hx1和hx2Hx2给出的滤波器输出。使用（16），（20），（30），（31）和（32），我们可以很容易地发现：

其中，Rx1x1、Rx2x2和Rx1x2分别是信号向量的x自相关和互相关。使用这些互谱估计器，我们现在能够获得Sα（w）使用第3.3节中描述的循环解调过程和图4所示的估计。程序如下：对于每个感兴趣的循环频率，我们移位x(n)一个步长±α/2以获得表达式（11）和（12），并使用这些频移版本的x(n)作为交叉PSD（33）、（34）或（35）中任一个的输入。其输出是在相应的Sa(w)循环频率a下的估计。

MATLAB代码：

%
%   File:      wosa_test01.m
​
​
%   Description:
%       Example of SCD estimation using WOSA for diferent filter sizes
%   Notes:
​
​
clear; clc;
​
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simulation Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N  = 16;        % Number of samples
fs = 1;         % Sampling Frequency
nsamp = 16;     % Pulse shaping oversampling
SNR = 10;       % SNR in DB
​
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Derived Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ts = 1/fs;
fd = 1/(8*ts);
td = 1/fd;
fc = 1/(5*ts);
snr = 10^(SNR/10);          % SNR as ratio
​
Nup = N*nsamp;      % Number of samples of the upsampled signal
​
t = [0:Nup-1]*ts;     % time vector
Nfft = pow2(nextpow2(Nup));
​
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BPSK Signal generation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sn = sign(rand(1,N)-0.5);   % Create BPSK signal
sn = rectpulse(sn,nsamp);
​
win = rectwin(Nfft)';
​
%% Carrier modulate
sn = sn.*cos(2*pi*fc*t);
​
%% Add some Gaussian noise
Eb = sum(abs(sn).^2)/length(sn);
No = Eb/snr;
noise_var = sqrt(No/2);
​
noise = noise_var * randn(1, length(sn));
​
sn_noisy = sn + noise;
​
L = length(sn_noisy);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Obtain the SCDs via different methods
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[S,f,a] = scd(sn_noisy,L/16,'
WOSA
',fs);
figure(2);
subplot(221);
contour(a,f,abs(S));
title('
BPSK Cyclic Spectrum (WOSA N/
16
)
','
FontSize
', 16);
xlabel('
\alpha
','
FontSize
', 16);
ylabel('
freq (Hz)
','
FontSize
', 16);
zlabel('
PSD
','
FontSize
', 16);
figure(5);
subplot(221);
mesh(a,f,abs(S));
title('
BPSK Cyclic Spectrum (WOSA N/
16
)
','
FontSize
', 16);
xlabel('
\alpha
','
FontSize
', 16);
ylabel('
freq (Hz)
','
FontSize
', 16);
zlabel('
PSD
','
FontSize
', 16);
​
[S,f,a] = scd(sn_noisy,L/8,'
WOSA
',fs);
figure(2);
subplot(222);
contour(a,f,abs(S));
title('
BPSK Cyclic Spectrum (WOSA N/
8
)
','
FontSize
', 16);
xlabel('
\alpha
','
FontSize
', 16);
ylabel('
freq (Hz)
','
FontSize
', 16);
zlabel('
PSD
','
FontSize
', 16);
figure(5);
subplot(222);
mesh(a,f,abs(S));
title('
BPSK Cyclic Spectrum (WOSA N/
8
)
','
FontSize
', 16);
xlabel('
\alpha
','
FontSize
', 16);
ylabel('
freq (Hz)
','
FontSize
', 16);
zlabel('
PSD
','
FontSize
', 16);
​
[S,f,a] = scd(sn_noisy,L/4,'
WOSA
',fs);
figure(2);
subplot(223);
contour(a,f,abs(S));
title('
BPSK Cyclic Spectrum (WOSA N/
4
)
','
FontSize
', 16);
xlabel('
\alpha
','
FontSize
', 16);
ylabel('
freq (Hz)
','
FontSize
', 16);
zlabel('
PSD
','
FontSize
', 16);
figure(5);
subplot(223);
mesh(a,f,abs(S));
title('
BPSK Cyclic Spectrum (WOSA N/
4
)
','
FontSize
', 16);
xlabel('
\alpha
','
FontSize
', 16);
ylabel('
freq (Hz)
','
FontSize
', 16);
zlabel('
PSD
','
FontSize
', 16);
​
[S,f,a] = scd(sn_noisy,L/2,'
WOSA
',fs);
figure(2);
subplot(224);
contour(a,f,abs(S));
title('
BPSK Cyclic Spectrum (WOSA N/
2
)
','
FontSize
', 16);
xlabel('
\alpha
','
FontSize
', 16);
ylabel('
freq (Hz)
','
FontSize
', 16);
zlabel('
PSD
','
FontSize
', 16);
figure(5);
subplot(224);
mesh(a,f,abs(S));
title('
BPSK Cyclic Spectrum (WOSA N/
2
)
','
FontSize
', 16);
xlabel('
\alpha
','
FontSize
', 16);
ylabel('
freq (Hz)
','
FontSize
', 16);
zlabel('
PSD
','
FontSize
', 16);

%
%   File:      scd.m
​
​
%   Description:
%       Estimates the Spectral Correlation Density of a signal using the
%       signal modulates method.
%
​
​
function [SCD F A] = scd(x,M,METHOD,fs,res)
​
% x      -  Process for which the SCD is to be estimated
% M      -  For WLS method means the segment size
%           For Capon/APES means the filter length
% METHOD -  Cross spectrum method = [ WLS, MLM, NMLM, APES ] 
% fs     -  Sampling frequency
% res    -  cycle resolution
%
% SCD    - The bi-frequency SCD
% F      - The frequency axis in Hz
% A      - The cyclic frequency axis in Hz
​
if exist('METHOD','var') == 0
    METHOD = 'wosa';
else
    %% Try to match one of the available methods
    methods_list = {'wosa','mvdr','nmvdr','apes'};
    indx = strmatch(lower(METHOD),methods_list);
​
    if isempty(indx) | length(indx) > 1,
        error('Ambiguous or invalid method specified.');
    end
    METHOD = methods_list{indx};
end
​
if exist('fs','var') == 0
    fs = 1;
end
​
if exist('res','var') == 0
    res = length(x);
end
​
x = x(:);  % make sure the vector is column ordered
​
ts = 1/fs;
​
N = length(x);      % Length of the process
n = [0:N-1]'*ts;
A = [0:res/2]*fs/res;   % Obtain the cyclic frequency axis from 0 to N or 0 to fs
​
if strcmp(METHOD,'
wosa
')
    for i = 1:length(A)
        u = x .* exp(-j*pi*A(i)*n);
        v = x .* exp(j*pi*A(i)*n);
        [SCD(:,i) f] = wls(u,v,M);
    end
elseif strcmp(METHOD,'
mvdr
')
    for i = 1:length(A)
        u = x .* exp(-j*pi*A(i)*n);
        v = x .* exp(j*pi*A(i)*n);
        [SCD(:,i) f] = capon(u,v,M);
    end
elseif strcmp(METHOD,'
nmvdr
')
    for i = 1:length(A)
        u = x .* exp(-j*pi*A(i)*n);
        v = x .* exp(j*pi*A(i)*n);
        [SCD(:,i) f] = capon2(u,v,M);
    end
elseif strcmp(METHOD,'
apes
')
    for i = 1:length(A)
        u = x .* exp(-j*pi*A(i)*n);
        v = x .* exp(j*pi*A(i)*n);
        [SCD(:,i) f] = apes(u,v,M);
    end
else
    error(['
Unknown method 
' METHOD]); 
end
​
​
F = f*fs;   % Convert to Hz
​
%% We calculate the SCD for cyclic frequencies from 0 to fs/2
%% so now following the symetry property we obtain the SCD for
%% negative frequencies by mirroring.
​
SCD = [fliplr(conj(SCD(:,2:end))) SCD(:,1:end-1)];
A = [-fliplr(A(:,2:end)) A(:,1:end-1)];

仿真结果