解调信号，Hilbert是一种有效的方法，但它对信号由要求，需要满足Bedrosian Product Theorem，才不会有大的误差，而且Hilbert变换的计算过程中涉及加窗，也会给结果带来误差。最明显的效应是在非整周期采样情况下，Hilbert变换会带来边缘飞翼，边缘误差很大。
实际上，实现解调的数学方法是很多的，如能量算子、Mandelstam法、Shekel法、Prony法等，而其中的能量算子法简单使用，能有效克服Hilbert解调的误差。
能量算子法的定义如下
其离散算法主要有两种：DESA-1和DESA-2，它们的计算公式分别如下
DESA-1:
DESA-2:
其中的ψ为Teager算子
知道了这些，既可以编写程序实现瞬时频率与包络的计算了，当然，成熟的算法需要考虑各种细节，如边缘的延拓等，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({'原始信号','瞬时频率','包络'})
可以看到，包络和瞬时频率都能正确的提取出来，只是没有顾及边缘效应，边端效果较差，但也只是差在几个点上。
网上查代码，还看到一个三点对称差分能量算子法，其算法思想如下
Phid.m ：
DEO3S.m：
三点对称差分：
求解该差分使用了传递函数方法：
系统：
其传递函数表示为：
将
作为该系统的输入，输出即为三点对称差分后的能量算子分离结果，即DESA中的：
DESA.m ：
对其代码稍作整理后，封装成如下函数
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))
|
发表于 2011-4-18 22:18
提升卡
变色卡
显身卡
