小波变换 分离影像低频部分_具有可调Q因子的小波变换应用实践部分（附代码）...

上一篇介绍了

iKnowledge：具有可调Q因子的小波变换（TQWT）​zhuanlan.zhihu.com

这一部分接着介绍其实际应用。一般我们能获取得到的信号都是有限长离散信号。为了将具有可调Q因子的小波变换（TQWT ）应用于有限长离散信号。需要对连续部分进行改造，具体如下。

设定原始信号为x(n)，信号长度为N，变换后的信号为y(n)，信号长度为N0。

1. 当N0<N时候，定义如下低通滤波

其中，X(k)和Y(k)表示x(n)和y(n)对应的离散傅里叶变换。对应示意图如下图所示。

2. 当N0>N，对应的低通滤波器公式为

对应的示意图如图所示

3. 当N1<N时候，定义如下高通滤波

4. 当N1>N时候，定义如下高通滤波

下图表示有限长度信号对应下的带通滤波。

其中，N0和N1取值为最近的偶数

V0和V1对应公式如下

对应matlab代码如下

function [V0,V1]=afb(X,N0,N1)
%afb: analysis filter bank
%converts a vector X into two vectors:
% V0 of length N0, V1 of length N1

% Need:
%     N0, N1, both even integers
%     2<=N0<=length(x)
%     2<=N1<=length(x)
%     N0+N1>length(x)

% Example 
% X=rand(1,20);
% [V0,V1]=afb(X,16,12);
% Y=sfb(V0,V1,20);
% max(abs(X-Y))

X=X(:).';      %X is row vector
N=length(X);

P=(N-N1)/2;
T=(N0+N1-N)/2-1;
S=(N-N0)/2;

% transition-band function
v=(1:T)/(T+1)*pi;
trans=(1+cos(v)).*sqrt(2-cos(v))/2;

% low-pass subband
V0=zeros(1,N0);     
V0(0+1)=X(0+1);     %dc term
V0((1:P)+1)=X((1:P)+1);   %pass-band (pos freq)
V0(P+(1:T)+1)=X(P+(1:T)+1).*trans;   %trans-band (pos freq)
V0(N0/2+1)=0;     %Nyquist freq
V0(N0-P-(1:T)+1)=X(N-P-(1:T)+1).*trans;     % trans-band(neg freq)
V0(N0-(1:P)+1)=X(N-(1:P)+1);     %pass-band (neg freq)

%high-pass subband
V1=zeros(1,N1);
V1(0+1)=0;
V1((1:T)+1)=X(P+(1:T)+1).*trans(T:-1:1);
V1(T+(1:S)+1)=X(P+T+(1:S)+1);
if rem(N,2)==0
    V1(N1/2+1)=X(N/2+1);
end
V1(N1-T-(1:S)+1)=X(N-P-T-(1:S)+1);
V1(N1-(1:T)+1)=X(N-P-(1:T)+1).*trans(T:-1:1);

输出Y0和Y1为

对应matlab代码如下

function Y=sfb(V0,V1,N)
% sfb: synthesis filter bank

N0=length(V0);
N1=length(V1);

S=(N-N0)/2;
P=(N-N1)/2;
T=(N0+N1-N)/2-1;

% transition-band function
v=(1:T)/(T+1)*pi;
trans=(1+cos(v)).*sqrt(2-cos(v))/2;

%add 1 to indices because Matlab indexing starts at 1 (not 0)
% low-pass subband
Y0= zeros(1,N);
Y0(0+1)=V0(0+1);   %dc term
Y0((1:P)+1)=V0((1:P)+1);   %pass-band (pos freq)
Y0(P+(1:T)+1)=V0(P+(1:T)+1).*trans;  %trans-band (pos freq)
Y0(P+T+(1:S)+1)=0;    %stop-band (pos freq)
if rem(N,2)==0
    Y0(N/2+1)=0;   %Nyquist freq (if N even)
end
Y0(N-P-T-(1:S)+1)=0;   %stop-band (neg freq)
Y0(N-P-(1:T)+1)=V0(N0-P-(1:T)+1).*trans;   %trans-band (neg freq)
Y0(N-(1:P)+1)=V0(N0-(1:P)+1);   % pass-band (neg freq)

% high-pass subband
Y1=zeros(1,N);
Y1(0+1)=0;   %dc term
Y1((1:P)+1)=0;   %stop-band(pos freq)
Y1(P+(1:T)+1)=V1((1:T)+1).*trans(T:-1:1);   % trans-band (pos freq)
Y1(P+T+(1:S)+1)=V1(T+(1:S)+1);   %pass-band (pos freq)
if rem(N,2)==0
    Y1(N/2+1)=V1(N1/2+1);   %Nyquist freq (if N even)
end
Y1(N-P-T-(1:S)+1)=V1(N1-T-(1:S)+1);   %pass-band (neg freq)
Y1(N-P-(1:T)+1)=V1(N1-(1:T)+1).*trans(T:-1:1);   % trans-band (neg freq)
Y1(N-(1:P)+1)=0;   % stop-band (neg freq)

Y=Y0+Y1;

前文提到的uDWT和逆uDWT代码如下

function X=uDFT(x)
% unitary DFT

N=length(x);
X=fft(x)/sqrt(N);
end

function x=uDFTinv(X)
% inverse unitary DFT

N=length(X);
x=sqrt(N)*ifft(X);
end

那么TQWT将信号x转化为小波变换对应公式如下

对应代码如下

function w=tqwt(x,Q,r,J)
% w=tqwt(x,Q,r,J)
% Tunable Q-factor wavelet transform (TQWT)
% INPUT
%    x-input signal (even length)
%    Q-Q-factor
%    r-ovrsampling rate (redundancy)
%    J-number of levels
% OUTPUT
%    w-wavelet coefficients
% NOTES
% w{j} is subband j for j=1,...,J+1
% w{1} is the highest-frequency subband signal 
% w{J+1} is the low-pass subband signal
% % Example (verify perfect reconstruction)
% Q=4;r=3;J=3;
% N=200;
% x=rand(1,N)
% w=tqwt(x,Q,r,J);
% y=itqwt(w,Q,r,N);
% max(abs(x-y))

% Reference: 'Wavelet Transform with Tunable Q-Factor'
check_params(Q,r,J);%

beta=2/(Q+1);
alpha=1-beta/r;
N=length(x);

Jmax=floor(log(beta*N/8)/log(1/alpha));
if J>Jmax
    J=Jmax;
    fprintf('Note: too many levels specified. n')
    if Jmax>0
        fprintf('Reduce levels to %dn',Jmax);
    else
        fprintf('Increase signal length. n');
    end
    w=[];
    return
end

X=uDFT(x);%
w=cell(1,J+1);

% Jstages:
for j=1:J
    N0=2*round(alpha^j*N/2);
    N1=2*round(beta*alpha^(j-1)*N/2);
    [X,W]=afb(X,N0,N1);%
    w{j}=uDFTinv(W);%
end
w{J+1}=uDFTinv(X);
end

在此之前，需要检查参数，对应函数代码如下

function check_params(Q,r,J)
[st,i]=dbstack(1);

%--------Check Q-----

if ~isscalar(Q)
    error('Error in %s: Q must be scalar n', st.file)
end

if ~isnumeric(Q)
    error('Error in %s: Q must be numeric n', st.file)
end

if Q<1
    error('Error in %s: Q must be greater than or equal to 1.0 n', st.file)
end

如果想要从TQWT变换逆变换求取原始信号，则对应公式如下

对应代码如下

function y=itqwt(w,Q,r,N)

check_params(Q,r);
beta=2/(Q+1);
alpha=1-beta/r;
J=length(w)-1;
Y=uDFT(w{J+1});

for j=J:-1:1
    W=uDFT(w{j});
    M=2*round(alpha^(j-1)*N/2);
    Y=sfb(Y,W,M);%
end
y=uDFTinv(Y);

以上就是所有关于TQWT的matlab代码，给一个简单的小例子，验证原始信号经过TQWT和TQWT逆变换之后的信号是否和原始信号一致。

clear all;clc;close all;
load data1.mat;
Q=2;r=3;J=15;

x=emg1;N=size(emg1,2);
w=tqwt(x,Q,r,J);
y=itqwt(w,Q,r,N);
max(abs(x-y));
plot(x,'r');hold on;plot(y,'b--')

到如下结果，发现完全一致。

欢迎关注公众号iKnowledge，这里将向您分享机器人、信号处理、人生感悟等等。