小波变换理解及测试

连续小波变换

%%%%%%%%%%%%%%%%%连续小波变换示 示例%%%%%%%%%%%%%%%%%%
clear;
clc;
fs=1024; %采样频率
f1=100;
f2=200;
t=0:1/fs:1;

s = cos(2*pi*f1*t).*(t>=0.1 & t<0.3) + sin(2*pi*f2*t).*(t>0.7);
% s=sin(2*pi*f1*t)+sin(2*pi*f2*t); %两个不同频率正弦信号合成的仿真信号
%%%%%%%%%%%%%%%%%小波时频图绘制%%%%%%%%%%%%%%%%%%
wavename='cmor3-3';
totalscal=100; %尺度序列的长度，即scal的长度
wcf=centfrq(wavename); %小波的中心频率
cparam=2*wcf*totalscal; %为得到合适的尺度所求出的参数
a=totalscal:-1:1;
scal=cparam./a; %得到各个尺度，以使转换得到频率序列为等差序列

coefs=cwt(s,scal,wavename); %得到小波系数
f=scal2frq(scal,wavename,1/fs); %将尺度转换为频率

figure;
imagesc(t,f,abs(coefs)); %绘制色谱图
colorbar;
xlabel('时间 t/s');
ylabel('频率 f/Hz');
title('小波时频图');

figure;
mesh(t,f,abs(coefs)); 
axis tight;    
colorbar;
xlabel('时间 t/s');
ylabel('频率 f/Hz');     
title(['小波时频图','(',num2str(wavename),')']);

结果如下：

原始信号输入：

变换以后结果：

离散小波变换

x=0:0.01:10;
noisdopp =10*sin(10*x);
noisdopp = [zeros(1, 1000) noisdopp zeros(1, 1000)];
[A,D] = dwt(noisdopp,'sym1');


figure;
ax1 = subplot(3,1,1);plot(noisdopp);grid; ylabel('noisdopp');
ax2 = subplot(3,1,2);plot(A);grid; ylabel('A');
ax3 = subplot(3,1,3);plot(D);grid; ylabel('D');
%linkaxes([ax1,ax2,ax3],'x');

clear all;
% The current extension mode is zero-padding (see dwtmode).

% Load original one-dimensional signal. 
load sumsin; 
s = sumsin; 
% Perform decomposition at level 3 of s using db1. 
[c,l] = wavedec(s,3,'sym1');
% Using some plotting commands,
% the following figure is generated.
figure;
subplot(3,1,1);plot(s);ylabel('sumsin');
subplot(3,1,2);plot(c);ylabel('c');
subplot(3,1,3);plot(l);ylabel('l');

clc
clear all
close all

% 当前延拓模式是补零
oldmode = dwtmode('zpd');

ts = 0.001;
fs = 1/ts;
t=0:ts:1;
N = length(t);
x = sin(2*pi*10*t) + sin(2*pi*50*t) + sin(2*pi*100*t) + 0.1*randn(1, length(t));

figure
plot(t,x);
xlabel('时间 t/s');
ylabel('幅值 A');

% 一维小波分解
[c,l] = wavedec(x,5,'db3');

% 重构第1～5层逼近系数
a5 = wrcoef('a',c,l,'db3',5);
a4 = wrcoef('a',c,l,'db3',4);
a3 = wrcoef('a',c,l,'db3',3);
a2 = wrcoef('a',c,l,'db3',2);
a1 = wrcoef('a',c,l,'db3',1);
% 显示逼近系数
figure
subplot(5,1,1);plot(t, a5);ylabel('a5');
subplot(5,1,2);plot(t, a4);ylabel('a4');
subplot(5,1,3);plot(t, a3);ylabel('a3');
subplot(5,1,4);plot(t, a2);ylabel('a2');
subplot(5,1,5);plot(t, a1);ylabel('a1');
xlabel('时间 t/s');

% 重构第1～5层细节系数
d5 = wrcoef('d',c,l,'db3',5);
d4 = wrcoef('d',c,l,'db3',4);
d3 = wrcoef('d',c,l,'db3',3);
d2 = wrcoef('d',c,l,'db3',2);
d1 = wrcoef('d',c,l,'db3',1);
% 显示细节系数
figure
subplot(5,1,1);plot(t, d5);ylabel('d5');
subplot(5,1,2);plot(t, d4);ylabel('d4');
subplot(5,1,3);plot(t, d3);ylabel('d3');
subplot(5,1,4);plot(t, d2);ylabel('d2');
subplot(5,1,5);plot(t, d1);ylabel('d1');
xlabel('时间 t/s');

% 第1层细节信号的包络谱
yh = hilbert(d3);
aabs = abs(yh);                 % 包络的绝对值
aabs = aabs - mean(aabs);
aangle = unwrap(angle(yh));     % 包络的相位
af = diff(aangle)/2/pi;         % 包络的瞬时频率，差分代替微分计算
% NFFT = 2^nextpow2(N);
NFFT = 2^nextpow2(1024*4);      % 改善栅栏效应
f = fs*linspace(0,1,NFFT);
YH = fft(yh, NFFT)/N;           % Hilbert变换复信号的频谱
A = fft(aabs, NFFT)/N;          % 包络的频谱

figure
plot(f,abs(YH))
title('原始复信号的Hilbert谱')
xlabel('频率f (Hz)')
ylabel('|YH(f)|')

figure
subplot(221)
plot(t, aabs')
title('包络的绝对值')
legend('包络分析结果', '真实包络')
subplot(222)
plot(t, aangle)
title('调制信号的相位')
subplot(223)
plot(t(1:end-1), af*fs)
title('调制信号的瞬时频率')
subplot(224)
plot(f,abs(A))
title('包络的频谱')
xlabel('频率f (Hz)')
ylabel('|A(f)|')

% 恢复延拓模式
dwtmode(oldmode);