diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml new file mode 100644 index 000000000..68dbe39cd --- /dev/null +++ b/.github/workflows/CompatHelper.yml @@ -0,0 +1,24 @@ +name: CompatHelper + +on: + schedule: + - cron: '00 * * * *' + +jobs: + CompatHelper: + runs-on: ${{ matrix.os }} + strategy: + matrix: + julia-version: [1.2.0] + julia-arch: [x86] + os: [ubuntu-latest] + steps: + - uses: julia-actions/setup-julia@latest + with: + version: ${{ matrix.julia-version }} + - name: Pkg.add("CompatHelper") + run: julia -e 'using Pkg; Pkg.add("CompatHelper")' + - name: CompatHelper.main() + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + run: julia -e 'using CompatHelper; CompatHelper.main()' diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml new file mode 100644 index 000000000..150f0c4e9 --- /dev/null +++ b/.github/workflows/TagBot.yml @@ -0,0 +1,12 @@ +name: TagBot +on: + schedule: + - cron: 0 * * * * +jobs: + TagBot: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v1 + - uses: JuliaRegistries/TagBot@v1 + with: + token: ${{ secrets.GITHUB_TOKEN }} diff --git a/.gitignore b/.gitignore index 05461b3b2..fedd299ab 100644 --- a/.gitignore +++ b/.gitignore @@ -29,3 +29,5 @@ /usr /docs/build + +Manifest.toml diff --git a/.travis.yml b/.travis.yml index b909c3557..dfe843061 100644 --- a/.travis.yml +++ b/.travis.yml @@ -2,9 +2,11 @@ language: julia os: - linux julia: - - 0.7 - 1.0 - 1.1 + - 1.2 + - 1.3 + - 1 - nightly matrix: allow_failures: diff --git a/Project.toml b/Project.toml index c9bc89e0a..86c9410da 100644 --- a/Project.toml +++ b/Project.toml @@ -1,9 +1,8 @@ name = "DSP" uuid = "717857b8-e6f2-59f4-9121-6e50c889abd2" -version = "0.5.2" +version = "0.6.2" [deps] -Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" @@ -15,9 +14,13 @@ IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e" [compat] Compat = ">= 1.3.0" -Polynomials = ">= 0.1.0" -julia = "0.7,1" -IterTools = ">= 1.0.0" +FFTW = "1.1" +OffsetArrays = "0.11" +Polynomials = "0.6" +Reexport = "0.2" +SpecialFunctions = "0.8, 0.9" +julia = "1" +IterTools = "1" [extras] DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" diff --git a/REQUIRE b/REQUIRE deleted file mode 100644 index 152430b11..000000000 --- a/REQUIRE +++ /dev/null @@ -1,7 +0,0 @@ -julia 0.7 -Polynomials 0.1.0 -Reexport -SpecialFunctions -FFTW -Compat 1.3 -IterTools 1.0 diff --git a/docs/Project.toml b/docs/Project.toml index cc17d9b3a..65d072b4e 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,3 +1,6 @@ [deps] DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" + +[compat] +Documenter = "0.24" diff --git a/docs/make.jl b/docs/make.jl index d58ef4396..7a7ecb52e 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,8 +1,10 @@ using Documenter, DSP +DocMeta.setdocmeta!(DSP, :DocTestSetup, :(using DSP); recursive=true) + makedocs( modules = [DSP], - format = :html, + format = Documenter.HTML(), sitename = "DSP.jl", pages = Any[ "Contents" => "contents.md", diff --git a/docs/src/util.md b/docs/src/util.md index 2f41a8fc6..c3a75dc88 100644 --- a/docs/src/util.md +++ b/docs/src/util.md @@ -1,16 +1,15 @@ # `Util` - utility functions -```@meta -DocTestSetup = quote - using DSP -end -``` + +!!! note + As of version 0.6.1 of DSP.jl, `fftfreq` and `rfftfreq` have been moved from + DSP.jl to [AbstractFFTs.jl](https://github.com/JuliaMath/AbstractFFTs.jl) + version 0.5 and above. You can also access these functions through + [FFTW.jl](https://github.com/JuliaMath/FFTW.jl) version 1.1 and above. ```@docs unwrap unwrap! hilbert -fftfreq -rfftfreq nextfastfft pow2db amp2db @@ -25,7 +24,3 @@ shiftsignal! alignsignals alignsignals! ``` - -```@meta -DocTestSetup = nothing -``` diff --git a/src/DSP.jl b/src/DSP.jl index 372f33de6..a0ac2cdae 100644 --- a/src/DSP.jl +++ b/src/DSP.jl @@ -1,5 +1,3 @@ -__precompile__() - module DSP using FFTW diff --git a/src/Filters/filt.jl b/src/Filters/filt.jl index 82ada9c02..f032efb22 100644 --- a/src/Filters/filt.jl +++ b/src/Filters/filt.jl @@ -237,9 +237,9 @@ end # Zero phase digital filtering by processing data in forward and reverse direction function iir_filtfilt(b::AbstractVector, a::AbstractVector, x::AbstractArray) zi = filt_stepstate(b, a) - zitmp = copy(zi) pad_length = 3 * (max(length(a), length(b)) - 1) t = Base.promote_eltype(b, a, x) + zitmp = similar(zi, t) extrapolated = Vector{t}(undef, size(x, 1)+pad_length*2) out = similar(x, t) @@ -312,9 +312,9 @@ end # Zero phase digital filtering for second order sections function filtfilt(f::SecondOrderSections{T,G}, x::AbstractArray{S}) where {T,G,S} zi = filt_stepstate(f) - zitmp = similar(zi) pad_length = 6 * length(f.biquads) t = Base.promote_type(T, G, S) + zitmp = similar(zi, t) extrapolated = Vector{t}(undef, size(x, 1)+pad_length*2) out = similar(x, t) diff --git a/src/Filters/stream_filt.jl b/src/Filters/stream_filt.jl index bbddc79a8..617649a04 100644 --- a/src/Filters/stream_filt.jl +++ b/src/Filters/stream_filt.jl @@ -105,7 +105,9 @@ mutable struct FIRArbitrary{T} <: FIRKernel{T} hLen::Int end -function FIRArbitrary(h::Vector, rate::Real, Nϕ::Integer) +function FIRArbitrary(h::Vector, rate_in::Real, Nϕ_in::Integer) + rate = convert(Float64, rate_in) + Nϕ = convert(Int, Nϕ_in) dh = [diff(h); zero(eltype(h))] pfb = taps2pfb(h, Nϕ) dpfb = taps2pfb(dh, Nϕ) @@ -117,7 +119,20 @@ function FIRArbitrary(h::Vector, rate::Real, Nϕ::Integer) inputDeficit = 1 xIdx = 1 hLen = length(h) - FIRArbitrary(rate, pfb, dpfb, Nϕ, tapsPerϕ, ϕAccumulator, ϕIdx, α, Δ, inputDeficit, xIdx, hLen) + FIRArbitrary( + rate, + pfb, + dpfb, + Nϕ, + tapsPerϕ, + ϕAccumulator, + ϕIdx, + α, + Δ, + inputDeficit, + xIdx, + hLen + ) end diff --git a/src/deprecated.jl b/src/deprecated.jl index bea953778..c3ddc2aa5 100644 --- a/src/deprecated.jl +++ b/src/deprecated.jl @@ -4,4 +4,8 @@ Base.@deprecate_binding TFFilter PolynomialRatio Base.@deprecate_binding BiquadFilter Biquad Base.@deprecate_binding SOSFilter SecondOrderSections +Base.@deprecate Frequencies(nreal::Int, n::Int, multiplier::Float64) FFTW.Frequencies(nreal, n, multiplier) +Base.@deprecate fftfreq(n::Int, fs::Real=1) FFTW.fftfreq(n, fs) +Base.@deprecate rfftfreq(n::Int, fs::Real=1) FFTW.rfftfreq(n, fs) + @deprecate conv2 conv diff --git a/src/dspbase.jl b/src/dspbase.jl index 270fec428..e36b3519a 100644 --- a/src/dspbase.jl +++ b/src/dspbase.jl @@ -170,51 +170,92 @@ function deconv(b::StridedVector{T}, a::StridedVector{T}) where T filt(b, a, x) end + + """ _zeropad!(padded::AbstractVector, u::AbstractVector, + padded_axes = axes(padded), data_dest::Tuple = (1,), data_region = CartesianIndices(u)) -Place the portion of `u` specified by `data_region` into `padded` starting at -location `data_dest`, and set the rest of `padded` to zero. This will mutate -`padded`. +Place the portion of `u` specified by `data_region` into `padded`, starting at +location `data_dest`. Sets the rest of `padded` to zero. This will mutate +`padded`. `padded_axes` must correspond to the axes of `padded`. -`padded` must start at index 1, but `u` can have arbitrary axes. """ -@inline function _zeropad!(padded::AbstractVector, - u::AbstractVector, - data_dest::Tuple = (1,), # Tuple to be consistent with ND case - data_region = CartesianIndices(u)) +@inline function _zeropad!( + padded::AbstractVector, + u::AbstractVector, + padded_axes = axes(padded), + data_dest::Tuple = (first(padded_axes[1]),), + data_region = CartesianIndices(u), +) datasize = length(data_region) - start_i = data_dest[1] - # Use axes to accommodate arrays that do not start at index 1 data_first_i = first(data_region)[1] - copyto!(padded, start_i, u, data_first_i, datasize) - - padded[1:start_i - 1] .= 0 - padded[start_i + datasize:end] .= 0 + dest_first_i = data_dest[1] + copyto!(padded, dest_first_i, u, data_first_i, datasize) + padded[first(padded_axes[1]):dest_first_i - 1] .= 0 + padded[dest_first_i + datasize : end] .= 0 padded end -# padded must start at index 1, but u can have arbitrary offset -@inline function _zeropad!(padded::AbstractArray{<:Any, N}, - u::AbstractArray{<:Any, N}, - data_dest::Tuple = ntuple(::Integer -> 1, N), - data_region = CartesianIndices(u)) where N - # Copy the data to the beginning of the padded array +@inline function _zeropad!( + padded::AbstractArray, + u::AbstractArray, + padded_axes = axes(padded), + data_dest::Tuple = first.(padded_axes), + data_region = CartesianIndices(u), +) fill!(padded, zero(eltype(padded))) - pad_dest_axes = UnitRange.(data_dest, data_dest .+ size(data_region) .- 1) - copyto!(padded, CartesianIndices(pad_dest_axes), u, data_region) + dest_axes = UnitRange.(data_dest, data_dest .+ size(data_region) .- 1) + dest_region = CartesianIndices(dest_axes) + copyto!(padded, dest_region, u, data_region) padded end -# Make the padded buffer, and then place u into it. See _zeropad! -@inline _zeropad(u, padded_size, args...) = - _zeropad!(similar(u, padded_size), u, args...) + +""" + _zeropad(u, padded_size, [data_dest, data_region]) + +Creates and returns a new base-1 index array of size `padded_size`, with the +section of `u` specified by `data_region` copied into the region of the new + array as specified by `data_dest`. All other values will be initialized to + zero. + +If either `data_dest` or `data_region` is not specified, then the defaults +described in [`_zeropad!`](@ref) will be used. +""" +function _zeropad(u, padded_size, args...) + padded = similar(u, padded_size) + _zeropad!(padded, u, axes(padded), args...) +end + +function _zeropad_keep_offset(u, padded_size, u_axes, args...) + ax_starts = first.(u_axes) + new_axes = UnitRange.(ax_starts, ax_starts .+ padded_size .- 1) + _zeropad!(similar(u, new_axes), u, args...) +end + +function _zeropad_keep_offset( + u, padded_size, ::NTuple{<:Any, Base.OneTo{Int}}, args... +) + _zeropad(u, padded_size, args...) +end + +""" + _zeropad_keep_offset(u, padded_size, [data_dest, dat_region]) + +Like [`_zeropad`](@ref), but retains the first index of `u` when creating a new +array. +""" +function _zeropad_keep_offset(u, padded_size, args...) + _zeropad_keep_offset(u, padded_size, axes(u), args...) +end +>>>>>>> master """ Estimate the number of floating point multiplications per output sample for an @@ -232,21 +273,19 @@ function optimalfftfiltlength(nb, nx) # Step through possible nffts and find the nfft that minimizes complexity # Assumes that complexity is convex first_pow2 = ceil(Int, log2(nb)) - prev_pow2 = ceil(Int, log2(nfull)) + max_pow2 = ceil(Int, log2(nfull)) prev_complexity = os_fft_complexity(2 ^ first_pow2, nb) pow2 = first_pow2 + 1 - while pow2 <= prev_pow2 + while pow2 <= max_pow2 new_complexity = os_fft_complexity(2 ^ pow2, nb) new_complexity > prev_complexity && break prev_complexity = new_complexity pow2 += 1 end - nfft = pow2 > prev_pow2 ? 2 ^ prev_pow2 : 2 ^ (pow2 - 1) + nfft = pow2 > max_pow2 ? 2 ^ max_pow2 : 2 ^ (pow2 - 1) - # L is the number of usable samples produced by each block - L = nfft - nb + 1 - if L > nx - # If L > nx, better to find next fast power + if nfft > nfull + # If nfft > nfull, then it's better to find next fast power nfft = nextfastfft(nfull) end @@ -284,6 +323,7 @@ end Transform the smaller convolution input to frequency domain, and return it in a new array. However, the contents of `buff` may be modified. """ + @inline function os_filter_transform!(A::NTuple{<:Any, AbstractArray{<:Real}}, p) fA = p * A[1] for a in A[2:end] @@ -306,6 +346,7 @@ end @inline function os_filter_transform!(buff::Tuple{AbstractArray{<:Complex}}, p!) copy(p! * buff[1]) + end """ @@ -332,6 +373,7 @@ end p! * buff # p! operates in place on buff buff .*= filter_fd ip! * buff # ip! operates in place on buff + end # Used by `unsafe_conv_kern_os!` to handle blocks of input data that need to be padded. @@ -376,6 +418,7 @@ function unsafe_conv_kern_os_edge!( out_stop, save_blocksize, sout_deficit, # How many output samples are missing if nffts > sout + tdbuff_axes, ) where N # Iterate over all possible combinations of edge dimensions for a number of # edges @@ -389,7 +432,7 @@ function unsafe_conv_kern_os_edge!( # on an edge. # # First make all entries equal to the center blocks: - copyto!(perimeter_range, 1, center_block_ranges, 1, N) + @inbounds copyto!(perimeter_range, 1, center_block_ranges, 1, N) # For the dimensions chosen to be on an edge (edge_dims), get the # ranges of the blocks that would need to be padded (lie on an edge) @@ -408,7 +451,7 @@ function unsafe_conv_kern_os_edge!( # The center region for non-edge dimensions has been specified above, # so finish specifying the region of the perimeter for this edge # block - for (i, dim) in enumerate(edge_dims) + @inbounds for (i, dim) in enumerate(edge_dims) perimeter_range[dim] = perimeter_edge_ranges[i] end @@ -416,7 +459,7 @@ function unsafe_conv_kern_os_edge!( block_region = CartesianIndices( NTuple{N, UnitRange{Int}}(perimeter_range) ) - for block_pos in block_region + @inbounds for block_pos in block_region # Figure out which portion of the input data should be transformed block_idx = convert(NTuple{N, Int}, block_pos) @@ -439,7 +482,7 @@ function unsafe_conv_kern_os_edge!( # Convolve portion of input - _zeropad!(tdbuff, u, pad_before .+ 1, data_region) + _zeropad!(tdbuff, u, tdbuff_axes, pad_before .+ 1, data_region) os_conv_block!(tdbuff, fdbuff, filter_fd, p, ip) # Save convolved result to output @@ -464,21 +507,22 @@ function unsafe_conv_kern_os_edge!( end # Assumes u is larger than, or the same size as, v +# nfft should be greater than or equal to 2*sv-1 function unsafe_conv_kern_os!(out, - u::AbstractArray{<:Any, N}, - v, - A, - su, - sv, - sout, - nffts) where N + u::AbstractArray{<:Any, N}, + v, + A, + su, + sv, + sout, + nffts) where N u_start = first.(axes(u)) out_axes = axes(out) out_start = first.(out_axes) out_stop = last.(out_axes) ideal_save_blocksize = nffts .- sv .+ 1 - # Number of samples that are "missing" if the valid portion of the - # convolution is smaller than the output + # Number of samples that are "missing" if the output is smaller than the + # valid portion of the convolution sout_deficit = max.(0, ideal_save_blocksize .- sout) # Size of the valid portion of the convolution result save_blocksize = ideal_save_blocksize .- sout_deficit @@ -486,6 +530,7 @@ function unsafe_conv_kern_os!(out, # Pre-allocation tdbuff, fdbuff, p, ip = os_prepare_conv(u, nffts) + tdbuff_axes = axes(tdbuff) # Transform the smaller filter _zeropad!(tdbuff, v) @@ -494,14 +539,16 @@ function unsafe_conv_kern_os!(out, for a in A)...), p) filter_fd .*= 1 / prod(nffts) # Normalize once for brfft - last_full_blocks = fld.(su, save_blocksize) # block indices for center blocks, which need no padding - center_block_ranges = UnitRange.(2, last_full_blocks) + first_center_blocks = cld.(sv .- 1, save_blocksize) .+ 1 + last_center_blocks = fld.(su, save_blocksize) + center_block_ranges = UnitRange.(first_center_blocks, last_center_blocks) + # block index ranges for blocks that need to be padded # Corresponds to the leading and trailing side of a dimension, or if there # are no center blocks, corresponds to the whole dimension - edge_ranges = map(nblocks, last_full_blocks) do nblock, lastfull - lastfull > 1 ? [1:1, lastfull + 1 : nblock] : [1:nblock] + edge_ranges = map(nblocks, first_center_blocks, last_center_blocks) do nblock, firstfull, lastfull + lastfull > 1 ? [1:firstfull - 1, lastfull + 1 : nblock] : [1:nblock] end all_dims = 1:N # Buffer to store ranges of indices for a single region of the perimeter @@ -548,14 +595,15 @@ function unsafe_conv_kern_os!(out, out_start, out_stop, save_blocksize, - sout_deficit) # How many output samples are missing if nffts > sout + sout_deficit, + tdbuff_axes) # How many output samples are missing if nffts > sout end tdbuff_region = CartesianIndices(tdbuff) # Portion of buffer with valid result of convolution valid_buff_region = CartesianIndices(UnitRange.(sv, nffts)) # Iterate over block indices (not data indices) that do not need to be padded - for block_pos in CartesianIndices(center_block_ranges) + @inbounds for block_pos in CartesianIndices(center_block_ranges) # Calculate portion of data to transform block_idx = convert(NTuple{N, Int}, block_pos) @@ -615,8 +663,6 @@ function _conv_kern_fft!(out, A::NTuple{<:Any, AbstractArray{T, N}}, Apad[1], CartesianIndices(UnitRange.(1, outsize))) end - - function _conv_fft!(out, A::Tuple{<:AbstractArray, <:AbstractArray}, S, outsize) os_nffts = map(optimalfftfiltlength, S[2], S[1]) if any(os_nffts .< outsize) @@ -640,8 +686,6 @@ function _conv_fft!(out, A, S, outsize) _conv_kern_fft!(out, A, outsize, nffts) end end - - # For arrays with weird offsets function _conv_similar(u, outsize, axes...) out_offsets = .+([first.(ax) for ax in axes]...) @@ -667,7 +711,8 @@ function _conv_sz(A, S) outsize = .+(S...) .- (length(S) - 1) out = _conv_similar(A, outsize) _conv!(out, A, S, outsize) -end + + # May switch argument order """ @@ -682,14 +727,11 @@ function _conv(A::AbstractArray...) maxnd = max([ndims(a) for a in A]...) return conv([cat(a, dims=maxnd) for a in A]...) end -# TODO: not the most concise notation, make base julia issue? + _conv(A::AbstractArray{<:Number, N}...) where {N} = conv(promote(A...)...) _conv(A::AbstractArray{<:Integer}...) = round.(Int, conv([float(a) for a in A]...)) function _conv(A::AbstractArray{<:BLAS.BlasFloat, N}...) where N sizes = size.(A) - # TODO: either delete comment or bring sorting back - # order = sortperm([p for p in prod.(sizes)], rev=true) - # _conv_sz(A[order], sizes[order]) _conv_sz(A, sizes) end @@ -763,9 +805,9 @@ function xcorr( padmode = check_padmode_kwarg(padmode, su, sv) if padmode == :longest if su < sv - u = _zeropad(u, sv) + u = _zeropad_keep_offset(u, sv) elseif sv < su - v = _zeropad(v, su) + v = _zeropad_keep_offset(v, su) end conv(u, dsp_reverse(conj(v), axes(v))) elseif padmode == :none diff --git a/src/periodograms.jl b/src/periodograms.jl index f4dad9e13..b33d4718e 100644 --- a/src/periodograms.jl +++ b/src/periodograms.jl @@ -6,6 +6,7 @@ using ..Util, ..Windows export arraysplit, nextfastfft, periodogram, welch_pgram, mt_pgram, spectrogram, power, freq, stft using FFTW +import FFTW: Frequencies, fftfreq, rfftfreq ## ARRAY SPLITTER @@ -212,17 +213,15 @@ Spectrogram object. Returns a tuple of frequency bin centers for a given Periodogram2 object. - -See also: [`fftfreq`](@ref), [`rfftfreq`](@ref) """ freq(p::TFR) = p.freq freq(p::Periodogram2) = (p.freq1, p.freq2) FFTW.fftshift(p::Periodogram{T,F}) where {T,F<:Frequencies} = - Periodogram(p.freq.nreal == p.freq.n ? p.power : fftshift(p.power), fftshift(p.freq)) + Periodogram(p.freq.n_nonnegative == p.freq.n ? p.power : fftshift(p.power), fftshift(p.freq)) FFTW.fftshift(p::Periodogram{T,F}) where {T,F<:AbstractRange} = p # 2-d FFTW.fftshift(p::Periodogram2{T,F1,F2}) where {T,F1<:Frequencies,F2<:Frequencies} = - Periodogram2(p.freq1.nreal == p.freq1.n ? fftshift(p.power,2) : fftshift(p.power), fftshift(p.freq1), fftshift(p.freq2)) + Periodogram2(p.freq1.n_nonnegative == p.freq1.n ? fftshift(p.power,2) : fftshift(p.power), fftshift(p.freq1), fftshift(p.freq2)) FFTW.fftshift(p::Periodogram2{T,F1,F2}) where {T,F1<:AbstractRange,F2<:Frequencies} = Periodogram2(fftshift(p.power,2), p.freq1, fftshift(p.freq2)) FFTW.fftshift(p::Periodogram2{T,F1,F2}) where {T,F1<:AbstractRange,F2<:AbstractRange} = p @@ -442,7 +441,7 @@ struct Spectrogram{T,F<:Union{Frequencies,AbstractRange}} <: TFR{T} time::FloatRange{Float64} end FFTW.fftshift(p::Spectrogram{T,F}) where {T,F<:Frequencies} = - Spectrogram(p.freq.nreal == p.freq.n ? p.power : fftshift(p.power, 1), fftshift(p.freq), p.time) + Spectrogram(p.freq.n_nonnegative == p.freq.n ? p.power : fftshift(p.power, 1), fftshift(p.freq), p.time) FFTW.fftshift(p::Spectrogram{T,F}) where {T,F<:AbstractRange} = p """ diff --git a/src/util.jl b/src/util.jl index 9ccae64b0..1fffd6f01 100644 --- a/src/util.jl +++ b/src/util.jl @@ -23,7 +23,6 @@ export hilbert, rmsfft, meanfreq, unsafe_dot, - polyfit, shiftin!, finddelay, shiftsignal, @@ -108,53 +107,6 @@ fftabs2type(::Type{Complex{T}}) where {T<:FFTW.fftwReal} = T fftabs2type(::Type{T}) where {T<:FFTW.fftwReal} = T fftabs2type(::Type{T}) where {T<:Union{Real,Complex}} = Float64 -## FREQUENCY VECTOR - -struct Frequencies <: AbstractVector{Float64} - nreal::Int - n::Int - multiplier::Float64 -end - -unsafe_getindex(x::Frequencies, i::Int) = - (i-1+ifelse(i <= x.nreal, 0, -x.n))*x.multiplier -function Base.getindex(x::Frequencies, i::Int) - (i >= 1 && i <= x.n) || throw(BoundsError()) - unsafe_getindex(x, i) -end -if isdefined(Base, :iterate) - function Base.iterate(x::Frequencies, i::Int=1) - i > x.n ? nothing : (unsafe_getindex(x,i), i + 1) - end -else - Base.start(x::Frequencies) = 1 - Base.next(x::Frequencies, i::Int) = (unsafe_getindex(x, i), i+1) - Base.done(x::Frequencies, i::Int) = i > x.n -end -Base.size(x::Frequencies) = (x.n,) -Base.step(x::Frequencies) = x.multiplier - -""" - fftfreq(n, fs=1) - -Return discrete fourier transform sample frequencies. The returned -Frequencies object is an AbstractVector containing the frequency -bin centers at every sample point. `fs` is the sample rate of the -input signal. -""" -fftfreq(n::Int, fs::Real=1) = Frequencies(((n-1) >> 1)+1, n, fs/n) - -""" - rfftfreq(n, fs=1) - -Return discrete fourier transform sample frequencies for use with -`rfft`. The returned Frequencies object is an AbstractVector -containing the frequency bin centers at every sample point. `fs` -is the sample rate of the input signal. -""" -rfftfreq(n::Int, fs::Real=1) = Frequencies((n >> 1)+1, (n >> 1)+1, fs/n) -FFTW.fftshift(x::Frequencies) = (x.nreal-x.n:x.nreal-1)*x.multiplier - # Get next fast FFT size for a given signal length const FAST_FFT_SIZES = [2, 3, 5, 7] """ diff --git a/test/REQUIRE b/test/REQUIRE deleted file mode 100644 index 4d5b425cf..000000000 --- a/test/REQUIRE +++ /dev/null @@ -1 +0,0 @@ -OffsetArrays \ No newline at end of file diff --git a/test/dsp.jl b/test/dsp.jl index 4f9654cdd..dc54a09d6 100644 --- a/test/dsp.jl +++ b/test/dsp.jl @@ -36,7 +36,7 @@ end @testset "conv" begin - @test optimalfftfiltlength(3, 1) == 3 # Should be at least the first input + @test optimalfftfiltlength(1, 3) == 1 # Should be at least the first input @testset "conv-1D" begin # Convolution @@ -177,7 +177,7 @@ end for numdim in Ns for elt in eltypes for nsmall in regular_nsmall - nfft = optimalfftfiltlength(nlarge, nsmall) + nfft = optimalfftfiltlength(nsmall, nlarge) test_os(elt, nlarge, nsmall, Val{numdim}(), nfft) end end @@ -200,6 +200,7 @@ end test_os(Float64, 25, 4, Val{1}(), 16) end + @testset "N-arg-conv" begin u = [1, 2, 3, 2, 1] v = transpose([6, 7, 3, 2]) @@ -261,6 +262,7 @@ end [1, 2, 3, 4], [5, 6, 7, 8]) end + end @testset "xcorr" begin @@ -297,6 +299,9 @@ end off_a = OffsetVector(a, -1:1) off_b = OffsetVector(b, 0:1) @test xcorr(off_a, off_b, padmode = :none) == OffsetVector(exp, -2:1) + + # Issue #288 + @test xcorr(off_a, off_b, padmode = :longest) == OffsetVector(vcat(0, exp), -3:1) end @testset "deconv" begin diff --git a/test/filt.jl b/test/filt.jl index 1f237cb78..e04b360bc 100644 --- a/test/filt.jl +++ b/test/filt.jl @@ -152,6 +152,19 @@ end @test x2_matlab ≈ filtfilt(b, a, x) end +# Make sure above doesn't crash for real coeffs & complex data. +@testset "filtfilt 1D Complex" begin + x2_matlab = readdlm(joinpath(dirname(@__FILE__), "data", "filtfilt_output.txt"),'\t') + + b = [ 0.00327922, 0.01639608, 0.03279216, 0.03279216, 0.01639608, 0.00327922] + a = [ 1. , -2.47441617, 2.81100631, -1.70377224, 0.54443269, -0.07231567] + x = readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_x.txt"),'\t') + + y = x .+ 1im .* randn(size(x, 1)) + + @test x2_matlab ≈ real.(filtfilt(b, a, y)) +end + ####################################### # # Test 2d filtfilt against matlab results diff --git a/test/resample.jl b/test/resample.jl index 2a297518e..4faf601a5 100644 --- a/test/resample.jl +++ b/test/resample.jl @@ -59,6 +59,16 @@ end idxUpper = idxLower*2 yDelta = abs.(y[idxLower:idxUpper].-yy[idxLower:idxUpper]) @test all(map(delta -> abs(delta) < 0.005, yDelta)) + + # Test Float32 ratio (#302) + f32_ratio = convert(Float32, ratio) + f32_y = resample(x, f32_ratio) + ty = range(0, stop =cycles, length =yLen) + yy = sinpi.(2*ty) + idxLower = round(Int, yLen/3) + idxUpper = idxLower*2 + yDelta = abs.(f32_y[idxLower:idxUpper].-yy[idxLower:idxUpper]) + @test all(map(delta -> abs(delta) < 0.005, yDelta)) end @testset "resample_filter" begin diff --git a/test/runtests.jl b/test/runtests.jl index a1d31b976..9f7a91218 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,5 @@ using DSP, FFTW, Test +import DSP: Frequencies, fftfreq, rfftfreq using Random: seed! diff --git a/test/util.jl b/test/util.jl index 44884a468..73af99f72 100644 --- a/test/util.jl +++ b/test/util.jl @@ -50,27 +50,6 @@ using Statistics: mean end @testset "fft helpers" begin - ## FFTFREQ - @test fftfreq(1) ≈ [0.] - @test fftfreq(2) ≈ [0., -1/2] - @test fftfreq(2, 1/2) ≈ [0., -1/4] - @test fftfreq(3) ≈ [0., 1/3, -1/3] - @test fftfreq(3, 1/2) ≈ [0., 1/6, -1/6] - @test fftfreq(6) ≈ [0., 1/6, 1/3, -1/2, -1/3, -1/6] - @test fftfreq(7) ≈ [0., 1/7, 2/7, 3/7, -3/7, -2/7, -1/7] - - @test rfftfreq(1) ≈ [0.] - @test rfftfreq(2) ≈ [0., 1/2] - @test rfftfreq(2, 1/2) ≈ [0., 1/4] - @test rfftfreq(3) ≈ [0., 1/3] - @test rfftfreq(3, 1/2) ≈ [0., 1/6] - @test rfftfreq(6) ≈ [0., 1/6, 1/3, 1/2] - @test rfftfreq(7) ≈ [0., 1/7, 2/7, 3/7] - - for n = 1:7 - @test fftshift(fftfreq(n)) ≈ fftshift([fftfreq(n);]) - end - @test meanfreq(sin.(2*π*10*(0:1e-3:10*π)),1e3) ≈ 10.0 rtol=1e-3 # nextfastfft