小波分析公式(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

推荐阅读

给力项目线报网会员可免费下载 加入会员
友情提醒: 请尽量登录购买,防止付款了不发货!
QQ交流群:226333560 站长微信:qgzmt2
温馨提示:本站提供的一切软件、教程和内容信息都来自网络收集整理,仅限用于学习和研究目的;不得将上述内容用于商业或者非法用途,否则,一切后果请用户自负,版权争议与本站无关。用户必须在下载后的24个小时之内,从您的电脑或手机中彻底删除上述内容。如果您喜欢该程序和内容,请支持正版,购买注册,得到更好的正版服务。我们非常重视版权问题,如有侵权请邮件与我们联系处理。敬请谅解!

给TA打赏
共{{data.count}}人
人已打赏
行业资讯

图片管理的软件(图片管理神器ACDSee究竟是款怎样的软件?)

2024-3-26 23:17:03

冒泡网资源行业资讯

抖音24小时无人直播音乐,不违规,不封号纯撸音浪,小白实操当天日入1000+

2024-3-26 23:33:24

0 条回复 A文章作者 M管理员
    暂无讨论,说说你的看法吧
个人中心
购物车
优惠劵
今日签到
有新私信 私信列表
搜索