Python y Julia

from scipy import signal
import matplotlib.pyplot as plt
import numpy as np
from scipy.fftpack import fft
from scipy.fftpack import fftshift

0.3s

using SySTools

2.2s

Generar una señal de prueba, un seno de 2Vrms cuya frecuencia se modula lentamente alrededor de 3kHz, contaminada con ruido blanco de magnitud decreciente exponencialmente, muestreado a 10kHz.

fs = 10e3
N  = 1e5
amp = 2 * np.sqrt(2)
noise_power = 0.01 * fs / 2
time = np.arange(N) / float(fs)
mod = 1000 * np.cos(2*np.pi * 0.25 * time) # moduladora
carrier = amp * np.sin(2*np.pi * 3e3 * time + mod)
noise = np.random.normal(scale=np.sqrt(noise_power), size=time.shape)
noise *= np.exp(-time / 5)
x = carrier + noise

0.0s

fs = 10e3
N  = 1e5
amp = 2 * sqrt(2)
noise_power = 0.01 * fs / 2
time = (0:N-1) ./ fs
mod = 1000 .* cos.(2pi * 0.25 .* time) # moduladora
carrier = amp .* sin.(2pi * 3e3 .* time .+ mod)
noise = randn(size(time)) .* sqrt(noise_power)
noise .*= exp.(-time ./ 5)
x = carrier .+ noise
;

1.2s

Compute and plot the spectrogram

window = signal.tukey(256) # ventana de Tukey de 256 muestras
f, t, Sxx = signal.spectrogram(x, fs, window)
fig, ax = plt.subplots()
ax.pcolormesh(t,f,Sxx)
ax.set(ylabel='Frequency [Hz]', xlabel='Time [sec]')
fig

10.7s

window = tukey(256, 0.5)
specplot(x, fs=fs, window=window, xaxis="Frequency (Hz)", yaxis="Time (s)")

47.5s

Ploteo en tiempo y frecuencia de la ventana utilizada

fig, ax = plt.subplots()
ax.plot(window)
ax.set(
  title  = "Tukey window",
	ylabel = "Amplitude",
	xlabel = "Sample",
	ylim   = [0, 1.1]
)
fig

0.6s

let
  ts = range(0, length=length(window), step=1/fs)

  plot(window, 
    title="Tukey window", 
    xaxis="Normalized magnitude (dB)",
    yaxis="Normalized frequency (cycles per sample)")
end

3.8s

A = fft(window,2048)/(len(window)/2.0)  # fft de 2048 puntos
freq = np.linspace(-0.5, 0.5, len(A))
response = 20 * np.log10(np.abs(fftshift(A/abs(A).max())))
fig, ax = plt.subplots()
ax.plot(freq,response)
ax.set_xlim(-0.5, 0.5)
ax.set_ylim(-120, 0)
ax.set(
  title  = "Frequency response of the Tukey window",
  ylabel = "Normalized magnitude [dB]",
  xlabel = "Normalized frequency [cycles per sample]"
)
fig

1.2s

let
  n = 2048

  a = fft(pad(window, n))
  a = abs.(a)
  a ./= maximum(abs.(a)) # Normalizo a 1
  
  response = 20 .* log10.(fftshift(a))
  
  frecs = range(-0.5, 0.5, length=n)
 

  plot(frecs, response, 
    title="Frequency response of the Tukey window", xaxis="f (Hz)",
    ylim=[-120, 0])
end

3.8s