小波分析公式（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, noi我爱线报网sedData, 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(N我爱线报网oise 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