能量算子+信号解调

解调信号，Hilbert是一种有效的方法，但它对信号由要求，需要满足Bedrosian Product Theorem，才不会有大的误差，而且Hilbert变换的计算过程中涉及加窗，也会给结果带来误差。最明显的效应是在非整周期采样情况下，Hilbert变换会带来边缘飞翼，边缘误差很大。

实际上，实现解调的数学方法是很多的，如能量算子、Mandelstam法、Shekel法、Prony法等，而其中的能量算子法简单使用，能有效克服Hilbert解调的误差。

知道了这些，既可以编写程序实现瞬时频率与包络的计算了，当然，成熟的算法需要考虑各种细节，如边缘的延拓等，DESA-1的一个简单代码实现如下

% Discrete Energy Separation Algorithm 1(DESA-1)for estimation purpose on an AM signal.

clc

clear

close all

% 原始信号

k=0.8;

w=2*pi/256;

q=2*pi/12;

n = 1:512;

A = 1+k*cos(w*n);

x = A.*cos(q*n);

x = x(:); % 转化为列向量

% 能量算子计算瞬时频率与包络

y = circshift(x,-1)-x;

ey = y.^2-circshift(y,1).*circshift(y,-1); % Teager算子

ex = x.^2-circshift(x,1).*circshift(x,-1); % Teager算子

w = acos(1-(ey+circshift(ey,1))./(4*ex)); % 瞬时频率

a = sqrt(ex./(1-(1-(ey+circshift(ey,1))./(4*ex)).^2)); % 包络

figure

plot([x,w,a])

legend({'原始信号','瞬时频率','包络'})

可以看到，包络和瞬时频率都能正确的提取出来，只是没有顾及边缘效应，边端效果较差，但也只是差在几个点上。

网上查代码，还看到一个三点对称差分能量算子法，其算法思想如下

对其代码稍作整理后，封装成如下函数

function [A,f] = DESA(x)

% 三点对称差分能量算子求解包络

% 算法详细描述：

% 2Ψ[x(n)]

% |a(n)|=-------------------------------------

% sqrt( Ψ[ x(n+1)-x(n-1) ] )

phix = DEO3S(x);

phix1 = DEO3S(circshift(x',-1)'-circshift(x',1)');

U = 2*phix;

L = abs(sqrt(phix1));

A=abs(U./L); % 应去除分母为0的点

A(A==inf)=0;

pha = 1-phix1/2./phix;

pha(pha>1 & pha<-1) = 0;

f = 0.5*acos(pha);

end

function f=DEO3S(x)

% 使用传递函数法求解三点对称差分能量算子

% x：行向量

% H(z)=z(1+2*z^-1+z^-2)/4;

Px=Phid(x);

Ns=length(Px);

w=2*pi*(-Ns/2:Ns/2)/Ns;

w=[w(1:Ns/2),w(Ns/2+1:Ns)]; % 去0点

z=exp(1i*w);

Hz=z.*(1+2*z.^-1+z.^-2)/4;

Xz=fft(Px,Ns);

Xz=[Xz(Ns/2+1:Ns),Xz(1:Ns/2)];%重新排列

Yw=Hz.*Xz;

Yw=[Yw(Ns/2+1:Ns),Yw(1:Ns/2)];

f=real(ifft(Yw));

end

function f=Phid(x)

% Teager能量算子：y=x(n)^2-x(n-1)*x(n+1)

x=x';

f=x.^2-circshift(x,1).*circshift(x,-1);

f=f';

end

对此函数的测试如下

% function testPDESA

clc

clear

close all

fs=5120;

Ns=2048;

t=(0:Ns-1)/fs;

x2=exp(-5*t).*cos(2*pi*20*t); % 实际包络

x=x2.*sin(2*pi*500*t);

[y,f]=DESA(x); % 能量算子解包络

yhilbert=abs(hilbert(x)); % Hilbert解包络

figure

subplot(2,1,1),

plot(t,[y; yhilbert; x])

title('解包络')

legend('能量算子包络','Hilbert包络','原始数据');

grid on

subplot(2,1,2),

plot(t,[abs(x2)-y;abs(x2)-yhilbert])

grid on

axis([0 0.4 -0.02 0.02]);

legend('能量算子包络误差','Hilbert包络误差');

title('解包络误差')

figure

plot(abs(f)/2/pi*fs)

title('能量算子解调得到的频率')

可以看到它求解得到的包络和Hilbert得到的包络相差无几，也可以看到能量算子法在一些过零点处有些误差，但还可以容忍，而Hilbert法的边端飞翼就确实太大了（能量算子只在边界一两个点上误差很大，毕竟能量算子中涉及离散差分，差分一次就少一点的信息，故造成边缘点的误差）；另外能量算子得到的瞬时频率貌似不照，其实如果用优化算法优化后，效果是可以得到很大提升的。

实际上，能量算子的实际应用可以利用一个matlab工具箱：hht_toolbox，它是一个希尔伯特黄变换的一个工具箱（不过貌似用的人不多，不如G.Rilling的代码好），里面有计算4个能量算子的函数：

desa：原始能量算子法

desa1：DESA-1法

desa1m：DESA-1法的平滑结果

desa2：DESA-2法

代码效果测试如下

% 测试各种能量算子算法的效果

clc

clear all

close all

ts = 0.001;

N = 200;

fs = 1/ts;

k = 0.8;

t = (0:N-1)*ts;

% 原始信号

fm = 10;

fc = 100;

Am = 1 + k*cos(2*pi*fm*t); % 包络

Am = 2+exp(-10*t).*cos(2*pi*fm*t); % 包络

Ac = cos(2*pi*fc*t);

y = Am.*Ac; % 信号调制

% 包络分析

yh = hilbert(y);

Ah = abs(yh); % 包络的绝对值

Ph = unwrap(angle(yh)); % 包络的相位

fh = diff(Ph)/2/pi; % 包络的瞬时频率，差分代替微分计算

% 能量算子解调

xn = y';

[W,A] = desa(xn, ts);

[W1,A1] = desa1(xn, ts);

[W1m,A1m] = desa1m(xn, ts);

[W2,A2] = desa2(xn, ts);

figure

subplot(211)

plot(t,[Ah',A,A1,A1m,A2,xn])

title('能量算子解调振幅')

legend({'Hilbert包络','desa包络','desa1包络','desa1m包络','desa2包络','原始信号'})

fh(end+1) = fh(end);

subplot(212)

plot(t,[fh'*fs,W,W1,W1m,W2])

title('能量算子解调频率')

legend({'Hilbert频率','desa频率','desa1频率','desa1m频率','desa2频率'})

% ylim([60,80])

figure

plot(t,[Ah'-Am',A-Am',A1-Am',A1m-Am',A2-Am'])

title('能量算子解调振幅的误差')

legend({'Hilbert包络','desa包络','desa1包络','desa1m包络','desa2包络'})

sum(sum((A1-Am').^2))

sum(sum((A1m-Am').^2))

sum(sum((A2-Am').^2))

可以看到，Hilbert仍然有很大的边缘飞翼效应，无论是包络还是解调频率；desa函数计算的包络和频率误差比较大；另三个函数的计算结果相差无几，且无边缘飞翼，误差也不大，如果非要做个好坏的比较的话，可以从解调的RMS值中看出

>> sum(sum((A1-Am').^2))

sum(sum((A1m-Am').^2))

sum(sum((A2-Am').^2))

ans =

0.0038

ans =

0.0039

ans =

0.0037

可以看到误差最小的是desa2，其次是desa1，最后是desa1-m，而之所以desa1-m要比desa1误差要大，是前者是对后的一个低通滤波，这相当于一个误差，但前者的效果上看会更平滑。

从上面可以看到能量算子的厉害之处，其实能量算子还可以再很多地方使用，具体可以参考文献《Signal processing using the Teager Energy Operator and other nonlinear operators》。

原文出处