【时频分析】基于原子模极小化的稀疏时频表示【附MATLAB代码】

作品简介

摘要

非平稳信号的分析和处理通常在时频域进行,而时频域是通过离散Gabor变换(DGT)得到的。由于加窗的影响,DGT得到的T-F表示会发生扩展,从而降低T-F域分析和处理的性能。为了获得良好的局部化的T-F表示,已经研究了使用N1-范数的稀疏感知方法。然而,他们需要将连续参数离散化到网格上,这会导致模型不匹配。在本文中,我们提出了一种方法来估计稀疏的T-F表示使用原子范数。原子范数使得稀疏优化无需连续参数的离散化。数值实验表明,该方法得到的T-F表示比传统方法更稀疏。

索引项-短时傅立叶变换(STFT),凸优化,原子范数,基追踪,半定规划。

1.引言

非平稳信号通常在时频域(T-F)中进行分析和处理。为了将时域信号转换到T-F域,通常使用短时傅立叶变换(STFT)或离散Gabor变换(DGT),这是由于其简单和易于理解的结构[1,2]。然而,由它获得的T-F表示由于分析信号的加窗而扩展。这种扩展可能会影响T-F域分析和处理的性能。

为了实现良好局部化的T-F表示,已经提出了许多方法[3-14]。重新分配和同步压缩方法旨在使用相位导数信息将扩展分量重新定位到原始位置[3-6]。它们的性能受到由于窗口化导致的组分混合的影响[15]。稀疏感知方法是强大的工具,可以抵抗这种分量和噪声的混合[7-14]。稀疏感知方法的目标是找到欠定系统的稀疏解。然而,典型的公式化的基础上,1-范数最小化涉及离散化的一个连续的参数到一个网格上。由于信号和预定义网格之间的模型失配,它可能会降低性能[16]。

最近,已经研究了使用原子范数的稀疏优化[17-20],并将其应用于许多应用,例如线谱估计[21,22],到达方向估计[23-25]和雷达中的目标定位[26]。原子范数不需要连续参数的离散化。因此,在稀疏T-F表示中引入原子范数可以得到更好的局部化T-F表示。在本文中,我们提出了一种估计方法的稀疏TF表示。在所提出的方法中,估计问题被制定为原子范数最小化的条件下,分析的时域信号可以重建。实验证实,所提出的方法提供了比传统方法更稀疏的T-F表示。

2.仿真结果

人工信号的T-F表示。每一列(从左到右)分别显示了通过DGT、重新分配方法、W1-范数最小化、窗口原子范数最小化和所提出的方法获得的T-F表示。下面一行在红色框中说明了这些放大图。

MATLAB部分代码

%% parameter
a = 2^4;
windowLen = 2^7;
bw = 0.04;
g = dpss(windowLen, bw*windowLen/2, 1);
​
%% artifical signal
L = 2^10;
t = (0:L-1)';
​
x1 = sin(2*pi*(90/L)*t);
x2 = chirp(t, 120/L, L-1, 330/L);
x3 = chirp(t, 250/L, L-1, 500/L,'quadratic');
X = [x1, x2, x3];
​
s = sum(X,2);
f = Padding(s, a);
​
%% Estimate T-F representation
rho = 0.01;
[T, X, nu] = sparseTFR_ANM(f, g, a, rho);
​
%% Display the estimated T-F representation
min_dB = 60;
figure('Position', [900 766 960 320])
​
subplot(1,2,1)
c = dgt(f, calcCanonicalDualWindow(g, a), a, 2^10);
imagesc([0, L-1], [0, 1], 20*log10(abs(c)./max(abs(c), [], 'all')))
axis xy
xlabel('Samples')
ylim([0, 0.5])
ylabel('Normalized frequency')
caxis([-min_dB, 0])
colorbar
title('DGT')
​
subplot(1,2,2)
display_ANMSpec(gca, T, X, a, min_dB, 1)
xlabel('Samples')
ylabel('Normalized frequency')
title('Proposed method')
​
set(gcf, 'InvertHardCopy', 'off')
% saveas(gcf,'../results/ANMSpec.png')
​
​
%% Padding function for periodicity
function y = Padding(x, shiftLen)
sigLen = length(x);
paddingLen = ceil(sigLen/shiftLen)*shiftLen - sigLen;
y = [x; zeros(paddingLen, 1)];
end
​
%% Discrete Gabor transform
function c = dgt(f, g, a, M)
N = length(f)/a;
ind = mod((0:N-1)*a + (0:length(g)-1)' - floor(length(g)/2), length(f)) + 1;
c = fft(f(ind).*g, M);
end

全部代码:

相关学习资料见面包多链接https://mbd.pub/o/author-a2mYl2tsbA==/work

欢迎加入我的知识星球:https://t.zsxq.com/11PZOV9jw,永久获取更多相关资料、代码。


创作时间: