fastlib / fCWT

The fast Continuous Wavelet Transform (fCWT) is a library for fast calculation of CWT.
Apache License 2.0
281 stars 54 forks source link

STFT Matlab code in the article #25

Open whip123 opened 1 year ago

whip123 commented 1 year ago

Hi, I am currently trying to replicate the comparison in the article. But I couldn't get the same result as the Benchmark results for synthetic data showed. Could you post the STFT code in the article here? Thanks a @fastlib @lschneiderbauer @felixdollack . (URGENT)

felixdollack commented 1 year ago

Hi @whip123,

I am not associated with this project. Just having a personal interest in it. You also posted commented here that you could not build the library for Matlab on Windows. In case you managed, could you comment in the other issue how you managed?

Best, Felix

fastlib commented 1 year ago

Dear @whip123,

The article in Nature Computational Science should have all the details about the parameters used for the STFT. Can you show me in which way you cannot match the benchmark results?

fastlib

whip123 commented 1 year ago

Dear @fastlib

For now, I couldn't normalize my result as shown in the article. Could you tell me how? Below is my Matlab code to replicate the results in the article

tic [s,f,t] = stft(signal,fs,Window=blackman(250,'period'),OverlapLength=200,FFTLength=512); STFTtime = toc sdb = mag2db(abs(s)); mesh(t,f,sdb); cc = max(sdb(:))+[-50 0]; ax = gca; ax.CLim = cc; view(2) ylim([0 120]) ylabel("Frequency (Hz)"); xlabel("Time (s)"); colorbar; image

Thank you. whip123

fastlib commented 1 year ago

We normalized every time-frequency matrix to a [0,1] range. As such, the power scale is transformed to an arbitary power unit that allows a fair comparison for the ridge extraction analysis we did subsequently. I hope this answers your question. For completeness I pasted the code below:

`function [res,ts,frqs] = stftwrapper(times,tstep,data,fs,f0,f1,tit,fftlength)

[s,freqs,t] = stft(data,fs,'Window',blackman(round(fs*0.5)),'OverlapLength',round(fs*0.40),'FFTLength',fftlength);

[v,s0]=min(abs(freqs-f0));
[v,s1]=min(abs(freqs-f1));

image(t,freqs(s0:s1,:),normalizeMatrix(abs(s(s0:s1,:))),"CDataMapping","direct")
res = normalizeMatrix(abs(s(s0:s1,:)));
set(gca, 'ydir', 'normal');
colormap jet;

ts = t;
frqs = freqs(s0:s1,:);

xlabel("Time (s)");
ylabel("Frequency (Hz)");

xticks(min(times+30):tstep:(max(times+30)))
xticklabels(min(times):tstep:(max(times)))

ttl = title(tit);
ttl.Units = 'Normalize'; 
ttl.FontSize = sp(6);
ttl.Position(1) = 0; % use negative values (ie, -0.1) to move further left
ttl.HorizontalAlignment = 'left';  
ttl.Color = 'white';

return`

whip123 commented 1 year ago

Hi @fastlib

Thanks for the code, could you give me all the parameters as well, such as fft length, f0, f1 and times.

Furthermore, I would also like to apologize for the mistakes in my previous post, as I found that I posted the wrong spectrogram after applying the code. Below shows the spectrogram I got after applying the code with 500ms Blackman window, and 400ms overlap. image However, in the spectrogram, there are 2 white lines, do u have any idea what I should change in the code? in order to get the same graph as the article shown

Thank you.