How to obtain time-averaged wavelet spectral density by the new version of CWT?

6 次查看（过去 30 天）
Shen Shen
Shen Shen 2024-1-28
评论： Shen Shen 2024-5-9
I know how to obtain time-averaged wavelet spectral density by the old version of CWT, and the result is similar to the power spectral density (PSD) by function PWELCH. But how to obtain wavelet spectral density by the new version of CWT? Given that the old version is no longger recommended. Thank you. Here is the code and signal I used.
load u.mat
Fs=800; % sampling frequency of signal 'u', in Hz;
% PSD by pwelch
SegOvlp=75;
windowL=1024;
window=hamming(windowL);
noverlap=fix(SegOvlp/100*windowL);
nfft=windowL;
[pu_F,f]=pwelch(u,window,noverlap,nfft,Fs);
Pu_F=pu_F*var(u)/trapz(f,pu_F); %normalization
% WSD by old version of CWT, which gives similar results to pwelch
wavename='gaus3';
wcf=centfrq(wavename);
fre=10.^(linspace(-0.11,2.6,100));
scal=wcf*Fs./fre;
coefs=cwt(u,scal,wavename);
pu_w=sum(abs(coefs).^2,2)/Fs;
Pu_w=pu_w*var(u)/abs(trapz(fre,pu_w)); %normalization
% figure
loglog(f,Pu_F,'k');
hold on
loglog(fre,Pu_w,'r');
xlabel('\itf/Hz');ylabel('\itPSD');
legend('PSD','WSD by old CWT')
set(gcf,'Units','centimeters','Position',[10 5 15 8]);
set(gcf,'Color','w');
set(gcf, 'PaperPositionMode', 'auto');

请先登录，再回答此问题。

采纳的回答

Wayne King
Wayne King 2024-5-6
Hi Shen, you can use cwtfilterbank and then use the timeSpectrum method (function). If you want scale-averaged power, then scaleSpectrum is also available.
Note both methods accept either the raw time series data, or the CWT coefficient matrix.
For example:
load kobe
% Construct filterbank object. Only required input is length of data
% but see the reference page for all the properties you can specify
fb = cwtfilterbank(SignalLength=length(kobe));
[tavgp,f] = fb.timeSpectrum(kobe);
semilogx(f,tavgp)
grid on
xlabel('Frequency (Hz)')
ylabel('Power')
% There is also a convenience plot syntax which may provide a convenient
% visualization for you. Call the methods with no output arguments
figure
fb.timeSpectrum(kobe);
% As I said, if you want you, you can compute the CWT coefficients separately and
% input those
cfs = fb.wt(kobe);
figure
fb.timeSpectrum(cfs);
I would recommend the following help pages:
cwtfilterbank
timeSpectrum
Note there are normalization options under timeSpectrum. Hope that helps,
Wayne
  1 个评论
Shen Shen
Shen Shen 2024-5-9
Ran in:
Wayne, thanks so much for your prompt reply. It does help. However, a new problem appears. And I will have to bother you.
I have tried 'timeSpectrum' with both my data and the 'kobe' data. I found that the spectral results obtained by 'timeSpectrum' and 'pwelch' are so different, whatever normalization options are used under timeSpectrum. The 'spectrum' by 'timeSpectrum' seems different from the power spectral density? I'm confused. Is it possible that you may give a clue on this?
load kobe
figure(1)
fb = cwtfilterbank(SignalLength=length(kobe));
fb.timeSpectrum(kobe);
figure(2)
pwelch(kobe);

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Continuous Wavelet Transforms 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by