Skip to content

Commit

Permalink
Add frequencystep method for narrowband spectra
Browse files Browse the repository at this point in the history
  • Loading branch information
dingraha committed Apr 12, 2024
1 parent 9b085d7 commit 934d0b3
Show file tree
Hide file tree
Showing 2 changed files with 29 additions and 3 deletions.
16 changes: 14 additions & 2 deletions src/narrowband.jl
Original file line number Diff line number Diff line change
Expand Up @@ -132,6 +132,18 @@ The frequencies are calculated using the `rfftfreq` function in the FFTW.jl pack
"""
@inline frequency(sm::AbstractNarrowbandSpectrum) = rfftfreq(inputlength(sm), samplerate(sm))


"""
frequencystep(sm::AbstractNarrowbandSpectrum)
Return the frequency step size `Δf` associated with the narrowband spectrum.
"""
@inline function frequencystep(sm::AbstractNarrowbandSpectrum)
m = inputlength(sm)
df = 1/(timestep(sm)*m)
return df
end

"""
PressureTimeHistory(sm::AbstractNarrowbandSpectrum, p=similar(halfcomplex(sm)))
Expand Down Expand Up @@ -459,7 +471,7 @@ end
@inline function Base.getindex(psa::PowerSpectralDensityAmplitude{false}, i::Int)
@boundscheck checkbounds(psa, i)
m = inputlength(psa)
df = 1/(timestep(psa)*m)
df = frequencystep(psa)
if i == 1
@inbounds hc_real = psa.hc[i]/m
return hc_real^2/df
Expand All @@ -473,7 +485,7 @@ end
@inline function Base.getindex(psa::PowerSpectralDensityAmplitude{true}, i::Int)
@boundscheck checkbounds(psa, i)
m = inputlength(psa)
df = 1/(timestep(psa)*m)
df = frequencystep(psa)
if i == 1 || i == length(psa)
@inbounds hc_real = psa.hc[i]/m
return hc_real^2/df
Expand Down
16 changes: 15 additions & 1 deletion test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@ using AcousticMetrics: p_ref
using AcousticMetrics: r2rfftfreq, rfft, rfft!, irfft, irfft!, RFFTCache, dft_r2hc, dft_hc2r
using AcousticMetrics: PressureTimeHistory
using AcousticMetrics: PressureSpectrumAmplitude, PressureSpectrumPhase, MSPSpectrumAmplitude, MSPSpectrumPhase, PowerSpectralDensityAmplitude, PowerSpectralDensityPhase
using AcousticMetrics: starttime, timestep, time, pressure, frequency, halfcomplex, OASPL
using AcousticMetrics: starttime, timestep, frequencystep, time, pressure, frequency, halfcomplex, OASPL
using AcousticMetrics: octave_fraction, band_start, band_end, cband_number
using AcousticMetrics: AbstractProportionalBands
using AcousticMetrics: ExactOctaveCenterBands, ExactOctaveLowerBands, ExactOctaveUpperBands
Expand Down Expand Up @@ -218,6 +218,8 @@ end
@test all(isapprox.(frequency(phase), freq_expected; atol=1e-12))
@test all(isapprox.(amp, amp_expected; atol=1e-12))
@test all(isapprox.(phase.*amp, phase_expected.*amp_expected; atol=1e-12))
@test frequencystep(amp) freq_expected[2]
@test frequencystep(phase) freq_expected[2]

# Make sure I can go from a PressureSpectrum to an PressureTimeHistory.
ap_from_ps = PressureTimeHistory(amp)
Expand Down Expand Up @@ -263,6 +265,8 @@ end
@test all(isapprox.(frequency(phase), freq_expected; atol=1e-12))
@test all(isapprox.(amp, amp_expected; atol=1e-12))
@test all(isapprox.(phase.*amp, phase_expected.*amp_expected; atol=1e-12))
@test frequencystep(amp) freq_expected[2]
@test frequencystep(phase) freq_expected[2]

# Make sure I can go from a PressureSpectrum to an PressureTimeHistory.
ap_from_ps = PressureTimeHistory(amp)
Expand Down Expand Up @@ -315,6 +319,8 @@ end
@test all(isapprox.(frequency(phase), freq_expected; atol=1e-12))
@test all(isapprox.(amp, amp_expected; atol=1e-12))
@test all(isapprox.(phase, phase_expected; atol=1e-12))
@test frequencystep(amp) freq_expected[2]
@test frequencystep(phase) freq_expected[2]

# Make sure I can go from a PressureSpectrum to an PressureTimeHistory.
ap_from_ps = PressureTimeHistory(amp)
Expand Down Expand Up @@ -363,6 +369,8 @@ end
@test all(isapprox.(frequency(phase), freq_expected; atol=1e-12))
@test all(isapprox.(amp, amp_expected; atol=1e-12))
@test all(isapprox.(phase.*amp, phase_expected.*amp_expected; atol=1e-12))
@test frequencystep(amp) freq_expected[2]
@test frequencystep(phase) freq_expected[2]

# Make sure I can convert a mean-squared pressure to a pressure spectrum.
psamp = PressureSpectrumAmplitude(amp)
Expand Down Expand Up @@ -430,6 +438,8 @@ end
@test all(isapprox.(frequency(phase), freq_expected; atol=1e-12))
@test all(isapprox.(amp, amp_expected; atol=1e-12))
@test all(isapprox.(phase, phase_expected; atol=1e-12))
@test frequencystep(amp) freq_expected[2]
@test frequencystep(phase) freq_expected[2]

# Make sure I can convert a mean-squared pressure to a pressure spectrum.
psamp = PressureSpectrumAmplitude(amp)
Expand Down Expand Up @@ -529,6 +539,8 @@ end
@test all(isapprox.(frequency(phase), freq_expected; atol=1e-12))
@test all(isapprox.(amp, amp_expected; atol=1e-12))
@test all(isapprox.(phase.*amp, phase_expected.*amp_expected; atol=1e-12))
@test frequencystep(amp) freq_expected[2]
@test frequencystep(phase) freq_expected[2]

# Make sure I can convert a PSD to a pressure spectrum.
psamp = PressureSpectrumAmplitude(amp)
Expand Down Expand Up @@ -597,6 +609,8 @@ end
@test all(isapprox.(frequency(phase), freq_expected; atol=1e-12))
@test all(isapprox.(amp, amp_expected; atol=1e-12))
@test all(isapprox.(phase, phase_expected; atol=1e-12))
@test frequencystep(amp) freq_expected[2]
@test frequencystep(phase) freq_expected[2]

# Make sure I can convert a PSD to a pressure spectrum.
psamp = PressureSpectrumAmplitude(amp)
Expand Down

0 comments on commit 934d0b3

Please sign in to comment.