Ricker小波及其频率切片小波变换

在前面的文章里曾经说过，小波基的选择根据领域的不同而不同，例如机械振动冲击信号分析中常用的morlet小波，结构损伤识别中常用的Marr小波，字典学习中效果较好的Laplace小波，图像处理中比较好用的双树复小波，还有地震信号处理中经常使用的Ricker小波，这篇文章简要描述下Ricker小波及其频率切片小波变换。

Ricker小波是地震信号分析处理中最常用的地震子波形式，其数学模型为：

下图为Ricker小波的示意图，主频是20 Hz，延迟时间为0.2 s，采样频率为200 Hz/s。可以看出，时域波形和频谱峰值分别对应信号时延和主频。此外，Ricker小波时域和频域都具有紧支性，且能量有限。对地震信号的谱分解的精确度越高，越能精确刻画地下介质层的精确构造，为地震资料解释、油气勘探等提供资料。

sampleLen = 1000;  % 波形采样点数
fs = 1000; % 采样频率/Hz
fm = 75;  % 中心频率/Hz
fig = 0;   % 是否画图
ap = 0.7;  % 子波峰值位置
no = -25;    % 加入白噪声/db，wgn函数的参数
timeDic = 1000*(0:1/fs:1/fs*(sampleLen-1)); % 转化为ms
% 生成Ricker小波及相应的加噪信号
[OriData, noisedData, noise]
waveData = OriData;
% 计算功率谱
N = length(waveData);%数据长度
%i = 0:N-1;
%t = i/fs;
f = (0:N-1)*fs/N;%横坐标频率
waveData2 = fft(waveData,N);%进行fft变换
mag = abs(waveData2);%幅值
Opower = mag.^2;
%含噪声数据计算
noisedData2 = fft(noisedData,N);%进行fft变换
mag = abs(noisedData2);%幅值
NDpower = mag.^2;
% 仅噪声数据计算
noise2 = fft(noise,N);%进行fft变换
mag = abs(noise2);%幅值
%f = (0:N-1)*fs/N;%横坐标频率
Npower = mag.^2;

绘制相应的波形及功率谱

h = figure;
titleStr = ['fs = ' num2str(fs) ',' 'fm=' num2str(fm) ',' 'no=' num2str(no), 'sampleLen = ' num2str(sampleLen) ];
subplot(3,2,1);
plot(timeDic,waveData);
title('Ricker wavelet');xlabel('Time /ms');ylabel('Amplitude /mV');
hold on;
box on;
set(gcf, 'color', 'w');
set(gca, 'FontName', 'Times New Roman');
xlim([1 N]);
subplot(3,2,2);
plot(f(1:length(f)/2), Opower(1:length(f)/2));
title('Ricker wavelet spectrum');
xlabel('Frequency /Hz');
ylabel('Amplitude');
hold on;
box on;xlim([1 fs/2]);
set(gcf, 'color', 'w');
set(gca, 'FontName', 'Times New Roman');

% 含噪声信号
subplot(3,2,3);
plot(timeDic,noisedData);
title('Noise Ricker wavelet');xlabel('Time /ms');ylabel('Amplitude /mV');
hold on;
box on;
set(gcf, 'color', 'w');
set(gca, 'FontName', 'Times New Roman');
xlim([1 N]);
subplot(3,2,4);
plot(f(1:length(f)/2), NDpower(1:length(f)/2));
title('Noise Ricker wavelet');
xlabel('Frequency /Hz');
ylabel('Amplitude');
hold on;
box on;xlim([1 fs/2]);
set(gcf, 'color', 'w');
set(gca, 'FontName', 'Times New Roman');

% 噪声
subplot(3,2,5);
plot(timeDic,noise);
title('Noise');xlabel('Time /ms');ylabel('Amplitude /mV');
hold on;
%grid on;
box on;
set(gcf, 'color', 'w');
set(gca, 'FontName', 'Times New Roman');
xlim([1 N]);
subplot(3,2,6);
plot(f(1:length(f)/2), Npower(1:length(f)/2));
title('Noise spectrum');
xlabel('Frequency /Hz');
ylabel('Amplitude');
box on;
xlim([1 fs/2]);
set(gcf, 'color', 'w');
set(gca, 'FontName', 'Times New Roman');
saveas(gcf, [titleStr '.fig']);
save([['Noised-Ricker-[' titleStr ']'] '.mat'], 'noisedData');
save([['Orignal-Ricker-[' titleStr ']'] '.mat'], 'OriData');

下面看一下Ricker小波的频率切片小波变换。频率切片小波汲取了短时傅里叶变换STFT 和连续小波变换CWT的优点，继承了小波函数的特性，可调的尺度系数使 频率切片小波变换时频分解的时频分辨率可控。而且，逆变换不再像小波变换那样依赖小波函数，因此，可以灵活地在时频空间上进行时频区域分割，通过重构分离出期望的信号分量。

中文文献有些还是描述不清，直接看英文原版的吧，建议多看几篇

s = waveData;
Fs = fs;
N=length(s);
s=s-sum(s)/N;%消去直流信号
f1=0.;% 低频界
f2=1000;%高频界
%[f1,f2] 是观测到的频率范围，可以更改
k1=fix(f1*N/Fs-0.5);
k2=fix(f2*N/Fs-0.5);
df=1;  %观测到的频率步长
if(k2>N/2+1) k2=N/2+1; end
fp= fix(k1:df:k2);   %%  fp if the observed frequency in discrete form

nl=length(fp);

kapa=sqrt(2)/2/0.025;  %kapa 是时频分辨率因子

Tn=512; %时域采样点数

[A] = GetFSWT(s,Fs,fp,kapa,Tn);   %频率切片小波变换

ss = GetInvFSWT(N,A,fp);%recover the signal from the band of fp;
B=sqrt(A.*conj(A));

mx= max(max(B));
B=fix(B*128/mx);
t= (0:Tn-1)*N/Fs/Tn;

figure;

set(gcf, 'PaperUnits', 'inches');
set(gcf, 'PaperPositionMode', 'manual');
lan = 2;
if lan==2
    aa = 18/(2.54*lan);
    bb = 12/(2.54*2);
else
    aa = 18/(2.54*lan);
    bb = 12/(2.54*2);
end
set(gcf, 'PaperPosition', [1 1 aa bb]);
set(gcf, 'color', 'w');

subplot(2,7,[1 5]);
plot((0:N-1),ss,'b');
ylabel('Amplitude (mV)', 'FontSize', 9, 'FontName', 'Arial');
box on;
axis tight;
xlim([0, 3000]);

subplot(2,7,[13 14]);
Y=fft(s,N);
z=Y.*conj(Y);
z=sqrt(z);
z(1)=0;
z(2)=0;
K1=fix(f2*N/Fs);
fp1=[0 : K1-1]/N*Fs;
plot(z(1:K1),fp1, 'b');
set(gca,'YDir','normal');
set(gca,'YAxisLocation','right');
ylabel('Frequency (Hz)', 'FontSize', 9, 'FontName', 'Arial');
xlabel('S', 'FontSize', 9, 'FontName', 'Arial');
box on;
axis tight;
colormap(jet);

subplot(2,7,[8 12]);
image(t*1000,fp*Fs/N,B');
set(gca,'ydir','normal');
xlabel('Time (ms)', 'FontSize', 9, 'FontName', 'Arial');
ylabel('Frequency (Hz)', 'FontSize', 9, 'FontName', 'Arial');
box on;
axis tight;
colormap(jet);

看下0到400Hz频率区间的切片变换

0到300Hz区间

看下带噪小波的频率切片小波变换

看下0到400Hz频率区间的切片变换

哈哈，还蛮好玩的，小波理论是真的优美。小波之美，美不胜收

最后看下白噪声信号的频率切片小波变换

详细代码见如下链接

https://mianbaoduo.com/o/bread/YpyblZhs