Plotting a spectrogram using DSP.jl

TheLateKronos · March 26, 2021, 1:03pm

I am looking to plot a spectrogram of a signal, but I am running into some issues. I have created a minimal working example, and will walk through it. What follows is a bit messy, but I think that I clear question emerges at the end, and that showing the messy-ness is a good way of showing my confusion, which is the problem:

So I load the packages, create a signal, and have a look to confirm that I have generated a signal that rises in frequency with time:

using DSP, Plots; plotlyjs()
fs = 44_100
ts = range(0, stop=5, step=1/fs)
signal = @. sin(2π*ts^2)
plot(ts, signal)

I did an FFT on it as well to make sure that things look allright:

fftplot(signal, fs) = plot( FFTW.fftfreq(length(signal), fs), abs.(fft(signal)./length(signal)), 
                            xguide="Frequency  / Hz", yguide="Magnitude"
)
fftplot(signal, fs)

So now I am looking to create a spectrogram. I do it, and plot it, as follows:

spec = spectrogram(signal, fs)
plot(spec.time, spec.freq, spec.power, xguide="Time  / s", yguide="Frequency  / Hz")

Can you see the skinny line down along the x-axis? That is the only content of the plot. So initially, the axis-scaling is all messed up. This has been the case for all 3 times I have tried, and every time, the data is squished up against the x-axis. But lets zoom in on the output we got:

Now, I can make out that we have a frequency that rises with time, but the plot does look quite bad. The y-ticks are gone, and I can not hover to see the y-values. I also would prefer it to be continuos in colour instead of controur-lines. I have tried calling the last plot-command with heatmap instead of plot, which produces the following:

I can zoom into this as well, but the resolution seems horrible, and I can not see what in m example should reduce :

In addition there is no colorbar, which I have not been able to add.

Soooooo after a while I realized that I have been specifying the n argument with my fs variable, which of course leads to issues. But when I sortet that out, I got the following:

spec = spectrogram(signal, 10; fs=fs)
heatmap(spec.time, spec.freq, spec.power, xguide="Time  / s", yguide="Frequency  / Hz")

The bars become shorter the larger my n-value is, which I don’t understand. This seems to be the reason why the default, which takes the signal length divided by 8 as n, ends up so “short”. Additionally, I no longer see the linear rise in frequency over time that I am expecting, so I think that this is even less correct somehow.

What is the correct way to plot a spectrogram produced by DPS.jl? Could a recipe for it be built in, so that new users do not have to run into the issues I have had?

PS: Switching to GR did the exact same ting, without the ability to zoom in.

TheLateKronos · March 26, 2021, 1:05pm

Something like what is produced in A really brief introduction to audio signal processing in Julia | seaandsailor would be great - unfortunatly, I don’t recognize the plotting commands, and don’t intent of swiching plotting library because of spectrograms.

rafael.guerra · March 26, 2021, 2:04pm

@LateKronos, two things to consider: (i) you are sampling at 44 kHz a sine wave that has 5 Hz maximum, that is why you need to zoom in; (ii) for better resolution on the spectrogram, you need to feed more parameters:

using DSP
using Plots; gr()

fs = 44_100  # Hz
ts = range(0, stop=5, step=1/fs)  # seconds
signal = @. sin(2π*1_000*ts^2)   # sweep over f = (1000*ts) Hz
#  plot(ts, signal)   # > 200K points, better to use InspectDR
n = length(signal)
nw = n÷50
spec = spectrogram(signal, nw, nw÷2; fs=fs)
heatmap(spec.time, spec.freq, spec.power, xguide="Time [s]", yguide="Frequency [Hz]")

TheLateKronos · March 26, 2021, 2:52pm

That helps so much! Playing around with the numbers of samples per fft of the spectrogram made things more clear. The frequency-axis also extends from 0 to fs/2, and so setting y-lims to what I could tell to be the interesting parts based on an FFT of the whole signal turned out to be very helpful as well. It took some tinkering, but the result is great

My code is now

fs = 44100
number_of_peices = 10_000
samples_per_fft = length(data)÷number_of_peices
spec = spectrogram(data, samples_per_fft; fs)
plot(spec.time, spec.freq, spec.power, ylims=(500, 5e3),
xguide="Time  / s", yguide="Frequency  / Hz")

, where data is a morse signal with some high frequency noise that fades in and out on top. The resulting spectrogram became

, and one can see the morse signl and the noise very clearly.

The only tiny issue now is that there is still noise where the spectrogram seems to say there is none, but perhaps it is just not loud enough. Thanks so much!

rafael.guerra · March 26, 2021, 8:38pm

NB:
The Wigner-Ville transform should achieve better time-frequency resolution. See this reference for a nice summary and comparison of different Time-Frequency transform methods.
For an old Julia implementation of the W-V transform, you may want to look here.

ericphanson · March 26, 2021, 9:13pm

You could also plot it in units of decibels (basically taking a log) by doing pow2db.(spec.power), which can often be a useful scaling for spectrograms (though I’m not sure it matters here).

baggepinnen · March 27, 2021, 4:15am

Look no further:

using SignalAnalysis, DSP
using SignalAnalysis.Units: kHz
fs = 44_100
ts = range(0, stop=5, step=1/fs)
s = @. sin(2π*ts^2)
aa = resample(s, 1000/fs)
plot(
    tfd(signal(aa, 1kHz), Wigner(nfft=4*1024, smooth=30, method=:CM1980, 
window=hamming)),
    yscale = :identity,
    ylims  = (0, 45),
    yticks = 0:2:45,
    title  = "Wigner-Ville distribution ",
) |> display

Notice that this is a heavy computation for large signals, hence the call to resample above. Before using the WVD, make sure you’ve studied and understand the caveats, like cross-term artifacts, associated with it.

Ashwani_Rathee · May 22, 2021, 11:34am

Just in case, if anybody want to plot the way its normally plotted

using SampledSignals, WAV
audio_test_url = "https://upload.wikimedia.org/wikipedia/commons/4/48/Piano-phrase.wav"
 y, fs = wavread("../Downloads/Piano-phrase.wav")
audio_test = SampleBuf(y,fs)
n = length(audio_test.data)
nw = n÷50
spec = spectrogram(mono(audio_test).data, nw, nw÷10; fs=fs)
heatmap(spec.time, spec.freq, pow2db.(spec.power), xguide="Time [s]", yguide="Frequency [Hz]")

Screenshot from 2021-05-22 17-03-44

theloni-monk · July 13, 2022, 12:31am

Hi, I am trying to recreate this but am getting an error that the spectrogram method you used is not defined. Did you include an additional package that is not listed in this snippet?

Ashwani_Rathee · July 13, 2022, 8:37am

Welcome to our community @theloni-monk !!
using DSP needs to be added

Thomas_Pilutti · March 31, 2023, 1:42pm

Hi - I’m working through the example @Ashwani_Rathee provided, and am getting error(s) with the spectrogram() command call. Using julia 1.8.3. Any advice?

julia> spec = spectrogram(mono(audio_test).data, nw, nw÷10; fs=fs)

ERROR: MethodError: no method matching spectrogram(::Matrix{Float64}, ::Int64, ::Int64; fs=44100.0f0)
Closest candidates are:
spectrogram(::AbstractVector{T}, ::Int64, ::Int64; onesided, nfft, fs, window) where T at ~/.julia/packages/DSP/wENjI/src/periodograms.jl:420
spectrogram(::AbstractVector{T}, ::Int64) where T at ~/.julia/packages/DSP/wENjI/src/periodograms.jl:420 got unsupported keyword argument “fs”
Stacktrace:
[1] top-level scope
@ REPL[37]:1

Daniel_Berge · March 31, 2023, 2:42pm

You are getting this error as spectrogram requires an AbstractVector. The mono function you are using from SampledSignals returns a [Nx1] matrix. You can fix this using the vec function returning that as a vector instead.

spec = spectrogram(vec(mono(audio_test).data), nw, nw÷10; fs=fs)

Thomas_Pilutti · March 31, 2023, 3:25pm

Ah, thanks @Daniel_Berge - I should have known that a nx1 matrix would not be interpreted as a vector. All good.