上一篇介绍了
iKnowledge:具有可调Q因子的小波变换(TQWT)zhuanlan.zhihu.com![0672441dd8bb4f5c063662b162283a12.png](https://i-blog.csdnimg.cn/blog_migrate/45bd7b56a14ff4d311dc8360d531c2bf.jpeg)
这一部分接着介绍其实际应用。一般我们能获取得到的信号都是有限长离散信号。为了将具有可调Q因子的小波变换(TQWT )应用于有限长离散信号。需要对连续部分进行改造,具体如下。
设定原始信号为x(n),信号长度为N,变换后的信号为y(n),信号长度为N0。
1. 当N0<N时候,定义如下低通滤波
![7eb52750064efb78fcc8c24a78605fe9.png](https://i-blog.csdnimg.cn/blog_migrate/892c6d2d0874bb8f7f4ca2a25793e49b.jpeg)
![0ef49ca390fc967a81c4972ec16b45c2.png](https://i-blog.csdnimg.cn/blog_migrate/ec83c58be908d6169609036ba06d7b9b.jpeg)
其中,X(k)和Y(k)表示x(n)和y(n)对应的离散傅里叶变换。对应示意图如下图所示。
![6105ed389d3a6f137c5427bfbf182d86.png](https://i-blog.csdnimg.cn/blog_migrate/658fab27702b3e45c63dbec819f6fc2d.jpeg)
2. 当N0>N,对应的低通滤波器公式为
![7eb52750064efb78fcc8c24a78605fe9.png](https://i-blog.csdnimg.cn/blog_migrate/892c6d2d0874bb8f7f4ca2a25793e49b.jpeg)
![0d23274dafafa25564e9b073ef776ff8.png](https://i-blog.csdnimg.cn/blog_migrate/c844aa135155f2a54cd3486b4e5dd0e0.jpeg)
对应的示意图如图所示
![8801a4b09f2d30f142fdaf9f57f43aea.png](https://i-blog.csdnimg.cn/blog_migrate/89797bc93d6493a13aae38e66100b9e6.jpeg)
3. 当N1<N时候,定义如下高通滤波
![7eb52750064efb78fcc8c24a78605fe9.png](https://i-blog.csdnimg.cn/blog_migrate/892c6d2d0874bb8f7f4ca2a25793e49b.jpeg)
![df5fbd770b6ca818deaa856ec01d29df.png](https://i-blog.csdnimg.cn/blog_migrate/2e69b138e52c9013efe52de5dcd64d6b.jpeg)
![8c952effad2dee1aa45c4410a23545d2.png](https://i-blog.csdnimg.cn/blog_migrate/77257282d7c85444d9eb40ef5ba30451.jpeg)
4. 当N1>N时候,定义如下高通滤波
![7eb52750064efb78fcc8c24a78605fe9.png](https://i-blog.csdnimg.cn/blog_migrate/892c6d2d0874bb8f7f4ca2a25793e49b.jpeg)
![656ddd4e5795e2d300b5f5327a214011.png](https://i-blog.csdnimg.cn/blog_migrate/694acbb45169fcf7523f48966ad0ccf3.jpeg)
![88dd31db0d131d6a152e2c4bee078da4.png](https://i-blog.csdnimg.cn/blog_migrate/90bc8d7ded23b59beb18f7b2aad43d35.jpeg)
下图表示有限长度信号对应下的带通滤波。
![7ff1f78b490ac40c95bec3915dc3f2cf.png](https://i-blog.csdnimg.cn/blog_migrate/b25306faaad2599407f7c114b7fd284d.jpeg)
其中,N0和N1取值为最近的偶数
![a69edfb87af65c0229eee8a3aae8efcc.png](https://i-blog.csdnimg.cn/blog_migrate/0e2e58e5feb8b1d1301379c444c79114.jpeg)
V0和V1对应公式如下
![f9d8139e038cd40f6b5eb38c9d2b565b.png](https://i-blog.csdnimg.cn/blog_migrate/4e2e96ebbcaa323ce4fb9989fa650ae7.jpeg)
![709df6f8f01aa6722766ee90950c02fe.png](https://i-blog.csdnimg.cn/blog_migrate/11e47994eef5ac363933fd53c040b4e3.jpeg)
对应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为
![dc11a5543bc248eb03b47ffb0536f92a.png](https://i-blog.csdnimg.cn/blog_migrate/91abcc63ad17609bc7874b88a4670ad8.jpeg)
![6c96d385505f07ab9cd4857d53f31562.png](https://i-blog.csdnimg.cn/blog_migrate/c3b7a3ed2859dea00cdf5db30ba5f9ee.jpeg)
对应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转化为小波变换对应公式如下
![a98f0c9e02b6c4f184b8c074c156f530.png](https://i-blog.csdnimg.cn/blog_migrate/df345904489c5533d0b186a6016c2ca2.jpeg)
对应代码如下
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变换逆变换求取原始信号,则对应公式如下
![547b60490651531e015931c9fb275cb7.png](https://i-blog.csdnimg.cn/blog_migrate/61fb22de216730e6985c0f213696f9a4.jpeg)
对应代码如下
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--')
到如下结果,发现完全一致。
![45c5af4f457bddde32338e7bf5f3a2e4.png](https://i-blog.csdnimg.cn/blog_migrate/b369dacb9c91a2a8dea5a375c512e461.jpeg)
欢迎关注公众号iKnowledge,这里将向您分享机器人、信号处理、人生感悟等等。
![3eab7882db032dbec4ebf37c519ccc39.png](https://i-blog.csdnimg.cn/blog_migrate/24603ef5cdbe3bf417c14708774c8bb6.jpeg)