Daubechies小波基的构造
% 此程序实现构造小波基
% periodic_wavelet.m
function ss=periodic_wavelet;
clear;clc;
% global MOMENT; % 消失矩阶数
% global LEFT_SCALET; % 尺度函数左支撑区间
% global RIGHT_SCALET; % 尺度函数右支撑区间
% global LEFT_BASIS; % 小波基函数左支撑区间
% global RIGHT_BASIS; % 小波基函数右支撑区间
% global MIN_STEP; % 最小离散步长
% global LEVEL; % 计算需要的层数(离散精度）
% global MAX_LEVEL; % 周期小波最大计算层数
[s2,h]=scale_integer;
[test,h]=scalet_stretch(s2,h);
wave_base=wavelet(test,h);
ss=periodic_waveletbasis(wave_base);
function [s2,h]=scale_integer;
% 本函数实现求解小波尺度函数离散整数点的值
% sacle_integer.m
MOMENT=10; % 消失矩阶数
LEFT_SCALET=0; % 尺度函数左支撑区间
RIGHT_SCALET=2*MOMENT-1; % 尺度函数右支撑区间
LEFT_BASIS=1-MOMENT; % 小波基函数左支撑区间
RIGHT_BASIS=MOMENT; % 小波基函数右支撑区间
MIN_STEP=1/512; % 最小离散步长
LEVEL=-log2(MIN_STEP); % 计算需要的层数(离散精度）
MAX_LEVEL=8; % 周期小波最大计算层数
h=wfilters('db10','r'); % 滤波器系数
h=h*sqrt(2); % FI(T)=SQRT(2)*SUM(H(N)*FI(2T-N)) N=0:2*MOMENT-1;
for i=LEFT_SCALET+1:RIGHT_SCALET-1
for j=LEFT_SCALET+1:RIGHT_SCALET-1
k=2*i-j+1;
if (k>=1&k a(i,j)=h(k); % 矩阵系数矩阵
else
a(i,j)=0;
end
end
end
[s,w]=eig(a); % 求特征向量，解的基
s1=s(:,1);
s2=[0;s1/sum(s1);0]; % 根据条件SUM(FI(T))=1,求解;
% 本函数实现尺度函数经伸缩后的离散值
% scalet_stretch.m
function [s2,h]=scalet_stretch(s2,h);
MOMENT=10; % 消失矩阶数
LEFT_SCALET=0; % 尺度函数左支撑区间
RIGHT_SCALET=2*MOMENT-1; % 尺度函数右支撑区间
LEFT_BASIS=1-MOMENT; % 小波基函数左支撑区间
RIGHT_BASIS=MOMENT; % 小波基函数右支撑区间
MIN_STEP=1/512; % 最小离散步长
LEVEL=-log2(MIN_STEP); % 计算需要的层数(离散精度）
MAX_LEVEL=8; % 周期小波最大计算层数
for j=1:LEVEL % 需要计算到尺度函数的层数
t=0;
for i=1:2:2*length(s2)-3 % 需要计算的离散点取值（0，1，2，3 -> 1/2, 3/2, 5/2)
t=t+1;
fi(t)=0;
for n=LEFT_SCALET:RIGHT_SCALET; % 低通滤波器冲击响应紧支撑判断
if ((i/2^(j-1)-n)>=LEFT_SCALET&(i/2^(j-1)-n) fi(t)=fi(t)+h(n+1)*s2(i-n*2^(j-1)+1); % 反复应用双尺度方程求解
end
end
end
clear s
n1=length(s2);
n2=length(fi);
for i=1:length(s2)+length(fi) % 变换后的矩阵长度
if (mod(i,2)==1)
s(i)=s2((i+1)/2); % 矩阵奇数下标为小波上一层(0,1,2,3)离散值
else
s(i)=fi(i/2); % 矩阵偶数下标为小波下一层(1/2,3/2,5/2)(经过伸缩变换后)的离散值
end
end
s2=s;
end
% 采用双尺度方程求解小波基函数 PSI(T)
% wavelet.m
function wave_base=wavelet(test,h);
MOMENT=10; % 消失矩阶数
LEFT_SCALET=0; % 尺度函数左支撑区间
RIGHT_SCALET=2*MOMENT-1; % 尺度函数右支撑区间
LEFT_BASIS=1-MOMENT; % 小波基函数左支撑区间
RIGHT_BASIS=MOMENT; % 小波基函数右支撑区间
MIN_STEP=1/512; % 最小离散步长
LEVEL=-log2(MIN_STEP); % 计算需要的层数(离散精度）
MAX_LEVEL=8; % 周期小波最大计算层数
i=0;
for t=LEFT_BASIS:MIN_STEP:RIGHT_BASIS; % 小波基支撑长度
s=0;
for n=1-RIGHT_SCALET:1-LEFT_SCALET % g(n)取值范围
if((2*t-n)>=LEFT_SCALET&(2*t-n) s=s+h(1-n+1)*(-1)^(n)*test((2*t-n)/MIN_STEP+1); % 计算任意精度的小波基函数值
end
end
i=i+1;
wave_base(i)=s;
end
采用多孔trous算法(undecimated wavelet transform)实现小波变换
clear;clc;
%% 1.生成信号
f=50; % 频率
fs=800; % 采样率
T=128; % 信号长度
n=1:T;
y=sin(2*pi*f*n/fs)+2*exp(-f*n/(4*fs)); % 信号
% y=circshift(y.',3).';
%% 2.正变换
l1=wfilters('db4','l')*sqrt(2)/2; % 参考低通滤波器
l1_zeros=[l1,zeros(1,T-length(l1))]; % 低通滤波器1
h1=wfilters('db4','h')*sqrt(2)/2; % 参考高通滤波器
h1_zeros=[h1,zeros(1,T-length(h1))]; % 高通滤波器1
low1=ifft(fft(y).*fft(l1_zeros)); % 低频分量1
high1=ifft(fft(y).*fft(h1_zeros)); % 高频分量1
l2=dyadup(l1); % 原滤波器插值
l2_zeros=[l2,zeros(1,T-length(l2))]; % 低通滤波器2
h2=dyadup(h1); % 原滤波器插值
h2_zeros=[h2,zeros(1,T-length(h2))]; % 高通滤波器2
low2=ifft(fft(low1).*fft(l2_zeros)); % 低频分量2
high2=ifft(fft(low1).*fft(h2_zeros)); % 高频分量2
%% 3.反变换
lr2=circshift(l2_zeros(end:-1:1).',1).'; % 重构低通滤波器2
hr2=circshift(h2_zeros(end:-1:1).',1).'; % 重构高通滤波器2
lr1=circshift(l1_zeros(end:-1:1).',1).'; % 重构低通滤波器1
hr1=circshift(h1_zeros(end:-1:1).',1).'; % 重构高通滤波器1
lowr=(ifft(fft(low2).*fft(lr2))+ifft(fft(high2).*fft(hr2))); % 重构低频分量1(lowr=low1)
r_s=(ifft(fft(lowr).*fft(lr1))+ifft(fft(high1).*fft(hr1))); % 重构源信号(r_s=y)
%% 4.绘图
figure(1);
plot(y);
title('源信号');
figure(2);
plot(low1,'r');
hold on;
plot(low2,'b');
legend('第一层低频','第二层低频');
figure(3);
plot(high1,'r');
hold on;
plot(high2,'b');
legend('第一层高频','第二层高频');
figure(4);
plot(low1,'r');
hold on;
plot(lowr,'b.');
legend('第一层低频','重构第一层低频');
figure(5);
plot(y,'r');
hold on;
plot(r_s,'b.');
legend('源信号','重构信号');
disp(norm(low1-lowr))
disp(norm(y-r_s))
平移变换平移法(cycle_spinning)消除gibbs效应
clear;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 1.原始信号
f=50; % 信号频率
fs=800; % 采样频率
N=128; % 采样点
% 信号赋值
n=1:N;
y=sin(2*pi*f*n/fs);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 2.噪声
noise=0.4*rand(1,128);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 3.染噪信号
y_noise=y+noise;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 4.硬消噪采用cycle_spinning技术
% 累加量
z5=zeros(1,N);
% 平移变换频移法
for i=1:N;
z=circshift(y_noise.',i-1).'; % 源信号右平移
[z1,z2]=lwt(z,'db3'); % 小波正变换
z2=zeros(1,N/2); % 高频分量全部为零（主要噪声，硬消噪）
z3=ilwt(z1,z2,'db3'); % 小波反变换
z4=circshift(z3.',-(i-1)).'; % 变换后信号左平移
z5=z5+z4/N; % 平均
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 5.显示
error=norm(y-z5)/norm(y); % 相对误差
figure(1);
subplot(2,1,1)
plot(y,'r');
legend('源信号');
subplot(2,1,2);
plot(y_noise);
legend('染噪信号');
figure(2);
subplot(2,1,1)
plot(y,'r');
legend('源信号');
title(error);
subplot(2,1,2);
plot(z5);
legend('消噪后信号');
