Signal processing

filtfilt(coef, x::SampledSignal)

Same as filtfilt, but retains sampling rate information.

resample(x::SampledSignal, rate[, coef])

Same as resample, but correctly handles sampling rate conversion.

filt(f, x::SampledSignal[, si])
filt(b, a, x::SampledSignal[, si])

Same as filt, but retains sampling rate information.

circconv(x, y)

Computes the circular convolution of x and y. Both vectors must be the same length.

circcorr(x)
circcorr(x, y)

Computes the circular correlation of x and y. Both vectors must be the same length.

compose(r, t, a; duration, fs)

Compose a signal from a reference signal and a list of arrival times and amplitudes.

Examples:

julia> x = cw(10kHz, 0.01, 44.1kHz)
julia> y1 = compose(x, [0.01, 0.03, 0.04], [1.0, 0.8, 0.6]; duration=0.05)
julia> y2 = compose(real(x), [10ms, 30ms, 40ms], [1.0, 0.8, 0.6]; duration=50ms)

decompose(r, x)
decompose(r, x, n; threshold, refine)

Decompose a signal as a sum of reference signal with some arrival times and amplitudes. If the number of signals n is non-zero, the function will limit the decomposition to the strongest n signals.

The resolution for arrival time is one sample. The threshold parameter controls the stopping criterion for the decomposition. If the relative change in the residual energy is less than threshold, the decomposition stops. The refine parameter controls whether the amplitudes are refined using an optimization algorithm. Without refinement, the amplitudes are simply the matched filter outputs at the arrival times.

Examples:

julia> x = compose(mseq(12), [0.1, 0.2], [1.0, 0.7]; duration=1.0, fs=8000)
julia> x += 0.1 * randn(size(x))
julia> decompose(mseq(12), x)
(time=[0.1, 0.2], amplitude=[1.00201, 0.70061], index=[801, 1601])

delay!(x, v)

Delay signal x by v units. The delay may be specified in samples or in time units (e.g. 2.3𝓈). The function supports delays of fractional or negative number of samples. The length of the returned signal is the same as the original, with any part of the signal that is shifted beyond the original length discarded.

Examples:

julia> x = signal([1.0, 2.0, 3.0], 10.0)
julia> delay!(x, 1)
SampledSignal @ 10.0 Hz, 3-element Vector{Float64}:
 0.0
 1.0
 2.0

julia> x = signal([1.0, 2.0, 3.0], 10.0)
julia> delay!(x, -1)
SampledSignal @ 10.0 Hz, 3-element Vector{Float64}:
 1.0
 2.0
 0.0

julia> x = signal([1.0, 2.0, 3.0, 2.0, 1.0], 10.0)
julia> delay!(x, 1.2)
SampledSignal @ 10.0 Hz, 5-element Vector{Float64}:
 -0.10771311593885118
  0.8061269251734184
  1.7488781412804602
  2.9434587943152617
  2.2755141401574472

julia> x = signal([1.0, 2.0, 3.0, 2.0, 1.0], 10.0)
julia> delay!(x, -0.11𝓈)
SampledSignal @ 10.0 Hz, 5-element Vector{Float64}:
  2.136038062052266
  2.98570052268015
  1.870314441196549
  0.9057002486847239
 -0.06203155191384361

delay(x, v)

Create a delayed version of signal x with v units of delay. See delay! for details.

demon(x; fs, downsample, method, cutoff)

Estimates DEMON spectrum. The output is highpass filtered with a cutoff frequency and downsampled. Supported downsampling methods are :rms (default), :mean and :fir.

downconvert(s, sps, fc)
downconvert(s, sps, fc, pulseshape; fs)

Converts passband signal centered around carrier frequency fc to baseband, and downsamples it by a factor of sps. If the pulseshape is specified to be nothing, downsampling is performed without filtering.

dzt(x, L, K)
dzt(x, L)

Compute the discrete Zak transform (DZT) of signal x with L delay bins and K Doppler bins. The length of signal x must be equal to LK. If K is not specified, it is assumed to be the length of x divided by L. Returns a K × L complex matrix.

If the frame rate of x if fs Sa/s, the delay bins are spaced at 1/fs seconds and the Doppler bins are spaced at fs/LK Hz. The DZT is scaled such that the energy of the signal is preserved, i.e., sum(abs2, x) ≈ sum(abs2, X).

For efficient computation of the DZT, K should product of small primes.

Examples:

julia> x = randn(ComplexF64, 4096)
4096-element Vector{ComplexF64}:
  :
julia> X = dzt(x, 64)
64×64 Matrix{ComplexF64}:
  :

findsignal(r, s)
findsignal(r, s, n; prominence, finetune, mingap, mfo)

Finds up to n strongest copies of reference signal r in signal s. The reference signal r should have a delta-like autocorrelation for this function to work well.

If the keyword parameter finetune is set to 0, approximate arrival times are computed based on a matched filter. If it is set to a positive value, an iterative optimization is performed to find more accruate arrival times within roughly finetune samples of the coarse arrival times.

The keyword parameter mingap controls the minimum number of samples between two arrivals.

Returns named tuple (time=t, amplitude=a, mfo=m) where t is a vector of arrival times, a is a vector of complex amplitudes of the arrivals, and m is the complex matched filter output (if keyword parameter mfo is set to true). The arrivals are sorted in ascending order of arrival times.

Examples:

julia> x = chirp(1000, 5000, 0.1, 40960; window=(tukey, 0.05))
julia> x4 = resample(x, 4)
julia> y4 = samerateas(x4, zeros(32768))
julia> y4[128:127+length(x4)] = real(x4)          # time 0.000775𝓈, index 32.75
julia> y4[254:253+length(x4)] += -0.8 * real(x4)  # time 0.001544𝓈, index 64.25
julia> y4[513:512+length(x4)] += 0.6 * real(x4)   # time 0.003125𝓈, index 129.0
julia> y = resample(y4, 1//4)
julia> y .+= 0.1 * randn(length(y))
julia> findsignal(x, y, 3)
(time = [0.000781, 0.001538, 0.003125], amplitude = [...], mfo=[])
julia> findsignal(x, y, 3; mfo=true)
(time = [0.000775, 0.001545, 0.003124], amplitude = [...], mfo=[...])

fir(n, f1)
fir(n, f1, f2; fs, method)

Designs a n-tap FIR filter with a passband from f1 to f2 using the specified method. If frame rate fs is not specified, f1 and f2 are given in normalized units (1.0 being Nyquist). If f1 is 0, the designed filter is a lowpass filter, and if f2 is nothing then it is a highpass filter.

This method is a convenience wrapper around DSP.digitalfilter.

Examples:

julia> lpf = fir(127, 0, 10kHz; fs=44.1kHz)   # design a lowpass filter
127-element Array{Float64,1}:
  ⋮

julia> hpf = fir(127, 10kHz; fs=44.1kHz)      # design a highpass filter
127-element Array{Float64,1}:
  ⋮

julia> bpf = fir(127, 1kHz, 5kHz; fs=44.1kHz) # design a bandpass filter
127-element Array{Float64,1}:
  ⋮

gmseq(m)
gmseq(m, θ)

Generates an generalized m-sequence of length 2^m-1 or tap specification m.

Generalized m-sequences are related to m-sequences but have an additional parameter θ. When θ = π/2, generalized m-sequences become normal m-sequences. When θ < π/2, generalized m-sequences contain a DC-component that leads to an exalted carrier after modulation. When θ is atan(√(2^m-1)), the m-sequence is considered to be period matched. Period matched m-sequences are complex sequences with perfect discrete periodic auto-correlation properties, i.e., all non-zero lag periodic auto-correlations are zero. The zero-lag autocorrelation is 2^m-1, where m is the shift register length.

This function currently supports shift register lengths between 2 and 30.

Examples:

julia> x = gmseq(3)         # generate period matched m-sequence
7-element Array{Complex{Float64},1}:
 0.3535533905932738 + 0.9354143466934853im
 0.3535533905932738 + 0.9354143466934853im
 0.3535533905932738 + 0.9354143466934853im
 0.3535533905932738 - 0.9354143466934853im
 0.3535533905932738 + 0.9354143466934853im
 0.3535533905932738 - 0.9354143466934853im
 0.3535533905932738 - 0.9354143466934853im

julia> x = gmseq(3, π/4)    # generate m-sequence with exalted carrier
7-element Array{Complex{Float64},1}:
 0.7071067811865476 + 0.7071067811865475im
 0.7071067811865476 + 0.7071067811865475im
 0.7071067811865476 + 0.7071067811865475im
 0.7071067811865476 - 0.7071067811865475im
 0.7071067811865476 + 0.7071067811865475im
 0.7071067811865476 - 0.7071067811865475im
 0.7071067811865476 - 0.7071067811865475im

goertzel(x, f, n; fs)

Detects frequency f in input signal using the Goertzel algorithm.

The detection metric returned by this function is the complex output of the Goertzel filter at the end of the input block. Typically, you would want to compare the magnitude of this output with a threshold to detect a frequency.

When a block size n is specified, the Goertzel algorithm in applied to blocks of data from the original time series.

hadamard(i, k)

Generate a vector with the entries of row i of a Walsh-Hadamard matrix of size 2ᵏ × 2ᵏ. Rows are numbered from 0 to 2ᵏ-1, so that i = 1 is the first non-trivial (not all ones) Hadamard sequence.

hadamard(k)

Generate a Walsh-Hadamard matrix of size 2ᵏ × 2ᵏ. Each row of the matrix is orthogonal to all other rows.

idzt(X)

Compute the inverse discrete Zak transform (DZT) of 2D K × L complex signal X with L delay bins and K Doppler bins. The length of the returned signal is LK.

See dzt for more details.

Examples:

julia> x = randn(ComplexF64, 4096)
4096-element Vector{ComplexF64}:
  :
julia> X = dzt(x, 64)
64×64 Matrix{ComplexF64}:
  :
julia> idzt(X) ≈ x
true

istft(Real, X; nfft, noverlap, window)

Compute the inverse short time Fourier transform (ISTFT) of one-sided STFT coefficients X which is based on segments with nfft samples with overlap of noverlap samples. Refer to DSP.Periodograms.spectrogram for description of the parameters.

For perfect reconstruction, the parameters nfft, noverlap and window in stft and istft have to be the same, and the windowing must obey the constraint of "nonzero overlap add" (NOLA). Implementation based on Zhivomirov 2019 and istft in scipy.

Examples:

julia> x = randn(1024)
1024-element Array{Float64,1}:
 -0.7903319156212055
 -0.564789077302601
  0.8621044972211616
  0.9351928359709288
  ⋮
  2.6158861993992533
  1.2980813993011973
 -0.010592954871694647

julia> X = stft(x, 64, 0)
33×31 Array{Complex{Float64},2}:
  ⋮

julia> x̂ = istft(Real, X; nfft=64, noverlap=0)
1024-element Array{Float64,1}:
 -0.7903319156212054
 -0.5647890773026012
  0.8621044972211612
  0.9351928359709288
  ⋮
  2.6158861993992537
  1.2980813993011973
 -0.010592954871694371

istft(Complex, X; nfft, noverlap, window)

Compute the inverse short time Fourier transform (ISTFT) of two-sided STFT coefficients X which is based on segments with nfft samples with overlap of noverlap samples. Refer to DSP.Periodograms.spectrogram for description of the parameters.

For perfect reconstruction, the parameters nfft, noverlap and window in stft and istft have to be the same, and the windowing must obey the constraint of "nonzero overlap add" (NOLA). Implementation based on Zhivomirov 2019 and istft in scipy.

Examples:

julia> x = randn(Complex{Float64}, 1024)
1024-element Array{Complex{Float64},1}:
  -0.5540372432417755 - 0.4286434695080883im
  -0.4759024596520576 - 0.5609424987802376im
                      ⋮
 -0.26493959584225923 - 0.28333817822701457im
  -0.5294529732365809 + 0.7345044619457456im

julia> X = stft(x, 64, 0)
64×16 Array{Complex{Float64},2}:
  ⋮

julia> x̂ = istft(Complex, X; nfft=64, noverlap=0)
1024-element Array{Complex{Float64},1}:
  -0.5540372432417755 - 0.4286434695080884im
 -0.47590245965205774 - 0.5609424987802374im
                      ⋮
  -0.2649395958422591 - 0.28333817822701474im
  -0.5294529732365809 + 0.7345044619457455im

mfilter(r, s)

Matched filter looking for reference signal r in signal s.

mseq(m)
mseq(m, θ)

Generates an m-sequence of length 2^m-1 or tap specification m.

m-sequences are sequences of +1/-1 values with near-perfect discrete periodic auto-correlation properties. All non-zero lag periodic auto-correlations are -1. The zero-lag autocorrelation is 2^m-1, where m is the shift register length.

This function currently supports shift register lengths between 2 and 30.

If specification m is provided, it should be a list of taps for the shift register. List of known m-sequence taps can be found in books or online.

Examples:

julia> x = mseq(3)                  # generate regular m-sequence
7-element Array{Float64,1}:
  1.0
  1.0
  1.0
 -1.0
  1.0
 -1.0
 -1.0

julia> x = mseq((1,3))              # generate m-sequence with specification (1,3)
7-element Array{Float64,1}:
  1.0
  1.0
  1.0
 -1.0
  1.0
 -1.0
 -1.0

pll(x)
pll(x, fc)
pll(x, fc, bandwidth; fs)

Phased-lock loop to track carrier frequency (approximately fc) in the input signal. If fc is not specified, the algorithm will attempt to track the dominant frequency.

rcosfir(β, sps)
rcosfir(β, sps, span)

Raised cosine filter.

removedc!(s; α)

DC removal filter. Parameter α controls the cutoff frequency. Implementation based on Lyons 2011 (3rd ed) real-time DC removal filter in Fig. 13-62(d).

See also: removedc

removedc(s; α)

DC removal filter. Parameter α controls the cutoff frequency. Implementation based on Lyons 2011 (3rd ed) real-time DC removal filter in Fig. 13-62(d).

See also: removedc!

rrcosfir(β, sps)
rrcosfir(β, sps, span)

Root-raised cosine filter.

upconvert(s, sps, fc)
upconvert(s, sps, fc, pulseshape; fs)

Converts baseband signal with sps symbols per passband sample to a real passband signal centered around carrier frequency fc.

whiten(x; nfft, noverlap, window, γ)

Spectral whitening of input signal x in the frequency domain. The parameters nfft, noverlap and window are required for the computation of STFT coefficients of x. Refer to DSP.Periodograms.spectrogram for description of the parameters. γ is a scaling or degree-of-flattening factor. The algorithm is based on Lee 1986.