Ufdjrw

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()