I am trying to get a third-octave frequency spectrum of a time signal.
The time signal is the acoustic pressure of rotational rotor noise which is harmonic. Its fundamental frequency is ff = n * N_b and for that reason, all frequencies should be multiples of ff.
Using fft I get the expected result: Multiples of the fundamental frequency are the relevant frequencies in the spectrum.
To get the third-octave frequency spectrum I wanted to use python acoustics, but the result of the function bandpass_third_octaves is not what I expected.
I expected the peaks from the fft frequency spectrum to be simply shifted to the third-octave centre frequencies with an adjusted amplitude. At least that is what I would like to get.
I think I interpret the output of bandpass_third_octaves incorrectly. Its output is a tuple containing the third-octave frequencies and a list of arrays which are supposed to contain the amplitude values as far as I know.
I currently use the arrays' maximum value as resulting amplitude because it works better than using the sum of it. It is possible that this interpretation is my mistake.
I would appreciate any help. There is no need for me to use python acoustics. Any solution to get the third-octave frequency spectrum would be awesome.
Edit: Using the mean instead of the maximum produces better results, but I am still not completely satisfied with it.
import matplotlib.pyplot as plt
import numpy as np
from scipy.sparse import spdiags
from scipy.signal import butter, lfilter, freqz, filtfilt, sosfilt
import acoustics.octave
#from acoustics.octave import REFERENCE
import acoustics.bands
from scipy.signal import hilbert
from acoustics.standards.iso_tr_25417_2007 import REFERENCE_PRESSURE
from acoustics.standards.iec_61672_1_2013 import (NOMINAL_OCTAVE_CENTER_FREQUENCIES,
NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES)
try:
from pyfftw.interfaces.numpy_fft import rfft
except ImportError:
from numpy.fft import rfft
def bandpass_filter(lowcut, highcut, fs, order=8, output='sos'):
"""Band-pass filter.
:param lowcut: Lower cut-off frequency
:param highcut: Upper cut-off frequency
:param fs: Sample frequency
:param order: Filter order
:param output: Output type. {'ba', 'zpk', 'sos'}. Default is 'sos'. See also :func:`scipy.signal.butter`.
:returns: Returned value depends on `output`.
A Butterworth filter is used.
.. seealso:: :func:`scipy.signal.butter`.
"""
nyq = 0.5 * fs
low = lowcut / nyq
high = highcut / nyq
output = butter(order / 2, [low, high], btype='band', output=output)
return output
def bandpass(signal, lowcut, highcut, fs, order=8, zero_phase=False):
"""Filter signal with band-pass filter.
:param signal: Signal
:param lowcut: Lower cut-off frequency
:param highcut: Upper cut-off frequency
:param fs: Sample frequency
:param order: Filter order
:param zero_phase: Prevent phase error by filtering in both directions (filtfilt)
A Butterworth filter is used. Filtering is done with second-order sections.
.. seealso:: :func:`bandpass_filter` for the filter that is used.
"""
sos = bandpass_filter(lowcut, highcut, fs, order, output='sos')
if zero_phase:
return _sosfiltfilt(sos, signal)
else:
return sosfilt(sos, signal)
class Frequencies:
"""
Object describing frequency bands.
"""
def __init__(self, center, lower, upper, bandwidth=None):
self.center = np.asarray(center)
"""
Center frequencies.
"""
self.lower = np.asarray(lower)
"""
Lower frequencies.
"""
self.upper = np.asarray(upper)
"""
Upper frequencies.
"""
self.bandwidth = np.asarray(bandwidth) if bandwidth is not None else np.asarray(self.upper) - np.asarray(
self.lower)
"""
Bandwidth.
"""
def __iter__(self):
for i in range(len(self.center)):
yield self[i]
def __len__(self):
return len(self.center)
def __str__(self):
return str(self.center)
def __repr__(self):
return "Frequencies({})".format(str(self.center))
def angular(self):
"""Angular center frequency in radians per second.
"""
return 2.0 * np.pi * self.center
class OctaveBand(Frequencies):
"""Fractional-octave band spectrum.
"""
def __init__(self, center=None, fstart=None, fstop=None, nbands=None, fraction=1,
reference=acoustics.octave.REFERENCE):
if center is not None:
try:
nbands = len(center)
except TypeError:
center = [center]
center = np.asarray(center)
indices = acoustics.octave.index_of_frequency(center, fraction=fraction, ref=reference)
elif fstart is not None and fstop is not None:
nstart = acoustics.octave.index_of_frequency(fstart, fraction=fraction, ref=reference)
nstop = acoustics.octave.index_of_frequency(fstop, fraction=fraction, ref=reference)
indices = np.arange(nstart, nstop + 1)
elif fstart is not None and nbands is not None:
nstart = acoustics.octave.index_of_frequency(fstart, fraction=fraction, ref=reference)
indices = np.arange(nstart, nstart + nbands)
elif fstop is not None and nbands is not None:
nstop = acoustics.octave.index_of_frequency(fstop, fraction=fraction, ref=reference)
indices = np.arange(nstop - nbands, nstop)
else:
raise ValueError("Insufficient parameters. Cannot determine fstart and/or fstop.")
center = acoustics.octave.exact_center_frequency(None, fraction=fraction, n=indices, ref=reference)
lower = acoustics.octave.lower_frequency(center, fraction=fraction)
upper = acoustics.octave.upper_frequency(center, fraction=fraction)
bandwidth = upper - lower
nominal = acoustics.octave.nominal_center_frequency(None, fraction, indices)
super(OctaveBand, self).__init__(center, lower, upper, bandwidth)
self.fraction = fraction
"""Fraction of fractional-octave filter.
"""
self.reference = reference
"""Reference center frequency.
"""
self.nominal = nominal
"""Nominal center frequencies.
"""
def __getitem__(self, key):
return type(self)(center=self.center[key], fraction=self.fraction, reference=self.reference)
def __repr__(self):
return "OctaveBand({})".format(str(self.center))
def bandpass_frequencies(x, fs, frequencies, order=8, purge=False, zero_phase=False):
""""Apply bandpass filters for frequencies
:param x: Instantaneous signal :math:`x(t)`.
:param fs: Sample frequency.
:param frequencies: Frequencies. Instance of :class:`Frequencies`.
:param order: Filter order.
:param purge: Discard bands of which the upper corner frequency is above the Nyquist frequency.
:param zero_phase: Prevent phase error by filtering in both directions (filtfilt)
:returns: Tuple. First element is an instance of :class:`OctaveBand`. The second element an array.
"""
if purge:
frequencies = frequencies[frequencies.upper < fs / 2.0]
return frequencies, np.array(
[bandpass(x, band.lower, band.upper, fs, order, zero_phase=zero_phase) for band in frequencies])
def bandpass_third_octaves(x, fs, frequencies=NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES, order=8, purge=False,
zero_phase=False):
"""Apply 1/3-octave bandpass filters.
:param x: Instantaneous signal :math:`x(t)`.
:param fs: Sample frequency.
:param frequencies: Frequencies.
:param order: Filter order.
:param purge: Discard bands of which the upper corner frequency is above the Nyquist frequency.
:param zero_phase: Prevent phase error by filtering in both directions (filtfilt)
:returns: Tuple. First element is an instance of :class:`OctaveBand`. The second element an array.
.. seealso:: :func:`octavepass`
"""
return bandpass_fractional_octaves(x, fs, frequencies, fraction=3, order=order, purge=purge, zero_phase=zero_phase)
def bandpass_fractional_octaves(x, fs, frequencies, fraction=None, order=8, purge=False, zero_phase=False):
"""Apply 1/N-octave bandpass filters.
:param x: Instantaneous signal :math:`x(t)`.
:param fs: Sample frequency.
:param frequencies: Frequencies. Either instance of :class:`OctaveBand`, or array along with fs.
:param order: Filter order.
:param purge: Discard bands of which the upper corner frequency is above the Nyquist frequency.
:param zero_phase: Prevent phase error by filtering in both directions (filtfilt)
:returns: Tuple. First element is an instance of :class:`OctaveBand`. The second element an array.
.. seealso:: :func:`octavepass`
"""
if not isinstance(frequencies, Frequencies):
frequencies = OctaveBand(center=frequencies, fraction=fraction)
return bandpass_frequencies(x, fs, frequencies, order=order, purge=purge, zero_phase=zero_phase)
def _sosfiltfilt(sos, x, axis=-1, padtype='odd', padlen=None, method='pad', irlen=None):
"""Filtfilt version using Second Order sections. Code is taken from scipy.signal.filtfilt and adapted to make it work with SOS.
Note that broadcasting does not work.
"""
from scipy.signal import sosfilt_zi
from scipy.signal._arraytools import odd_ext, axis_slice, axis_reverse
x = np.asarray(x)
if padlen is None:
edge = 0
else:
edge = padlen
# x's 'axis' dimension must be bigger than edge.
if x.shape[axis] <= edge:
raise ValueError("The length of the input vector x must be at least " "padlen, which is %d." % edge)
if padtype is not None and edge > 0:
# Make an extension of length `edge` at each
# end of the input array.
if padtype == 'even':
ext = even_ext(x, edge, axis=axis)
elif padtype == 'odd':
ext = odd_ext(x, edge, axis=axis)
else:
ext = const_ext(x, edge, axis=axis)
else:
ext = x
# Get the steady state of the filter's step response.
zi = sosfilt_zi(sos)
# Reshape zi and create x0 so that zi*x0 broadcasts
# to the correct value for the 'zi' keyword argument
# to lfilter.
#zi_shape = [1] * x.ndim
#zi_shape[axis] = zi.size
#zi = np.reshape(zi, zi_shape)
x0 = axis_slice(ext, stop=1, axis=axis)
# Forward filter.
(y, zf) = sosfilt(sos, ext, axis=axis, zi=zi * x0)
# Backward filter.
# Create y0 so zi*y0 broadcasts appropriately.
y0 = axis_slice(y, start=-1, axis=axis)
(y, zf) = sosfilt(sos, axis_reverse(y, axis=axis), axis=axis, zi=zi * y0)
# Reverse y.
y = axis_reverse(y, axis=axis)
if edge > 0:
# Slice the actual signal from the extended signal.
y = axis_slice(y, start=edge, stop=-edge, axis=axis)
return y
rho = 1.2
a = 340
N_b = 1
R = 1
r_H = 10
A = np.pi*R**2
TA = 287
M_H = 0.3
w = M_H*a/R
n = w/(2*np.pi)
t = np.linspace(0,0.8,num=40000)
az = t*2*np.pi*n*N_b
sin = np.sin(az)
cos = np.cos(az)
#Thickness Noise
F_H = R/r_H
F_E = 0.00012875807653441588 #Bestimmt für den Propeller aus Paper
T1 = ((3-M_H*sin)*sin)/((1-M_H*sin)**3)
T2 = (M_H*(cos**2))/(10*(1-M_H*sin)**4)
T3 = 50 + 39*(M_H**2) - 45*M_H*sin - 11*(M_H**2)*(sin**2) + 12* (M_H**3) *sin - 18*(M_H**3)*(sin**3)
T_M = ((M_H**3)/12)*(-T1 + T2 * T3)
p_T = 0.5 * rho * a**2 * F_H * F_E * T_M
#Loading Noise
F_T = (TA/ (rho * a**2))**(3/2) * (1 / (60 * np.sqrt(2) * N_b))
L = 60 + 30 * M_H**2 * cos**2 - 120 * M_H * sin - 30 * M_H**3 * sin * cos**2 + 80 * M_H**2 * sin**2 + 9 * M_H**4 * sin**2 * cos**2 - 20 * M_H**3 * sin**3
L_M = cos * (1 - M_H * sin)**(-3) * L
p_L = 0.5 * rho * a**2 * F_H * F_T * L_M
#Total
p_total = p_T + p_L
plt.figure(1)
plt.plot(t, p_total)
plt.title('Signal in time domain')
plt.xlabel('time [s]')
plt.ylabel('acoustic pressure [Pa]')
#fundamental frequency
ff = n*N_b
print('ff',ff)
#Sampling frequency
T = t[1] - t[0]
f_s = 1/T
print('fs',f_s)
#Trying to get the one third octave frequency spectrum
test = bandpass_third_octaves(p_total, f_s,frequencies=NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES,order=8,purge=False,zero_phase = True)
a_l = list()
i = 0
while i < 34:
a = max(test[1][i])
a_l.append(a)
i+=1
f = NOMINAL_THIRD_OCTAVE_CENTER_FREQUENCIES
plt.figure(2)
plt.bar(f, np.abs(a_l))
plt.title('Supposed one third octave spectrum of the time signal')
plt.xlabel('frequency [Hz]')
plt.ylabel('acoustic pressure [Pa]')
plt.xlim(0,100)
#FFT of the time signal p_total
N = p_total.size
f = np.linspace(0, 1/T, N)
f_scaled = f[:N // 2]
p_total -= np.mean(p_total)
fft = np.fft.fft(p_total)
fft_scaled = np.abs(fft)[:N // 2] * 2 / N
plt.figure(3)
plt.bar(f_scaled, fft_scaled)
plt.title('Signal in frequency domain')
plt.xlabel('frequency [Hz]')
plt.ylabel('acoustic pressure [Pa]')
plt.xlim(0,100)
plt.show()
