From 6233f1c71bcfe17fc40b053dd3a27048b5c2759f Mon Sep 17 00:00:00 2001 From: eric-brandao Date: Thu, 5 Oct 2023 14:16:15 -0300 Subject: [PATCH] python 3.11 compatibility, some plot updates, and deconv with zero padding --- pytta/_h5utils.py | 4 +- pytta/_plot.py | 13 ++++-- pytta/classes/signal.py | 100 +++++++++++++++++++++++++++------------- pytta/generate.py | 6 +-- 4 files changed, 81 insertions(+), 42 deletions(-) diff --git a/pytta/_h5utils.py b/pytta/_h5utils.py index f1dfb85..6af8b3d 100644 --- a/pytta/_h5utils.py +++ b/pytta/_h5utils.py @@ -31,14 +31,14 @@ def int_parser(attr): if isinstance(attr, (np.int16, np.int32, np.int64, - np.int, + int, np.int0)): return int(attr) else: return attr def float_parser(attr): - if isinstance(attr, (np.float, + if isinstance(attr, (float, np.float16, np.float32, np.float64)): diff --git a/pytta/_plot.py b/pytta/_plot.py index 1bee0f0..c0a0970 100644 --- a/pytta/_plot.py +++ b/pytta/_plot.py @@ -394,7 +394,7 @@ def freq(sigObjs, smooth, xLabel, yLabel, yLim, xLim, title, decimalSep): if xLabel is None: xLabel = 'Frequency [Hz]' if yLabel is None: - yLabel = 'Magnitude' + yLabel = 'Magnitude [dB]' if title is None: title = '' @@ -426,10 +426,6 @@ def freq(sigObjs, smooth, xLabel, yLabel, yLim, xLim, title, decimalSep): ax.set_ylabel(yLabel, fontsize=16) ax.legend(loc='best', fontsize=12) - if xLim is None: - xLim = [data['minFreq'], data['maxFreq']] - ax.set_xlim(xLim) - if yLim is None: yLim = [np.nanmin(yLims[:,0]), np.nanmax(yLims[:,1])] margin = (yLim[1] - yLim[0]) / 20 @@ -443,6 +439,13 @@ def freq(sigObjs, smooth, xLabel, yLabel, yLim, xLim, title, decimalSep): ax.xaxis.set_major_formatter(ticker.ScalarFormatter(useOffset=True)) for item in (ax.get_xticklabels() + ax.get_yticklabels()): item.set_fontsize(14) + freq_oct_list = [16, 31.5, 63, 125, 250, 500, 1000, 2000, 4000, 8000, 16000] + freq_oct_list_str = [str(item) for item in freq_oct_list] + ax.set_xticks(freq_oct_list) + ax.set_xticklabels(freq_oct_list_str) + if xLim is None: + xLim = [data['minFreq'], data['maxFreq']] + ax.set_xlim(xLim) return fig diff --git a/pytta/classes/signal.py b/pytta/classes/signal.py index c3ca83d..c3e869b 100644 --- a/pytta/classes/signal.py +++ b/pytta/classes/signal.py @@ -414,6 +414,18 @@ def crop(self, startTime, endTime): startIdx = np.where(self.timeVector >= startTime)[0][0] self.timeSignal = self.timeSignal[startIdx:endIdx,:] + def zero_padding(self, num_zeros = 1000): + """ zero-padding the signal object + """ + newTimeSignal = np.zeros((self.timeSignal.shape[0] + num_zeros, + self.timeSignal.shape[1])) + # fill initial part with original + newTimeSignal[:self.timeSignal.shape[0], + :self.timeSignal.shape[1]] = self.timeSignal + # Make the SignalObj equal to original + zeros + zp_signal = SignalObj(newTimeSignal, 'time', self.samplingRate) + return zp_signal + def mean(self): print('DEPRECATED! This method will be renamed to', ':method:``.channelMean()``', @@ -1182,6 +1194,10 @@ class ImpulsiveResponse(_base.PyTTaObj): * "linear": Computes using the spectral division of the signals; + + * "linear_zp": + zero pads input and output (double size) and + computes using the spectral division of the signals; * "H1": Uses power spectral density Ser divided by See, with @@ -1284,7 +1300,7 @@ def __init__(self, excitation=None, recording=None, "passing as parameter the 'excitation' " + "and 'recording' signals, or a calculated " + "'ir'.") - # Zero padding + # Initial zero padding elif excitation.numSamples > recording.numSamples: print("Zero padding on IR calculation!") excitSig = excitation.timeSignal @@ -1302,7 +1318,8 @@ def __init__(self, excitation=None, recording=None, np.zeros((recSig.shape[0], excitSig.shape[1])) newTimeSignal[:excitSig.shape[0], :excitSig.shape[1]] = \ excitSig - excitation.timeSignal = newTimeSignal + excitation.timeSignal = newTimeSignal + self._methodInfo = {'method': method, 'winType': winType, 'winSize': winSize, 'overlap': overlap} self._systemSignal = self._calculate_tf_ir(excitation, @@ -1396,6 +1413,40 @@ def _to_dict(self): out = {'methodInfo': self.methodInfo} return out + def _calculate_regu_spk(self, inputSignal, outputSignal): + """ Computes regularized spectrum + """ + data = _make_pk_spectra(inputSignal.freqSignal) + outputFreqSignal = _make_pk_spectra(outputSignal.freqSignal) + freqVector = inputSignal.freqVector + b = data * 0 + 10**(-200/20) # inside signal freq range + a = data * 0 + 1 # outside signal freq range + minFreq = np.max([inputSignal.freqMin, outputSignal.freqMin]) + maxFreq = np.min([inputSignal.freqMax, outputSignal.freqMax]) + # Calculate epsilon + eps = self._crossfade_spectruns(a, b, + [minFreq/np.sqrt(2), + minFreq], + freqVector) + if maxFreq < np.min([maxFreq*np.sqrt(2), + inputSignal.samplingRate/2]): + eps = self._crossfade_spectruns(eps, a, + [maxFreq, + maxFreq*np.sqrt(2)], + freqVector) + eps = \ + eps \ + * float(np.max(np.abs(outputFreqSignal)))**2 \ + * 1/2 + C = np.conj(data) / \ + (np.conj(data)*data + eps) + C = _make_rms_spectra(C) + C = SignalObj(C, + 'freq', + inputSignal.samplingRate, + signalType='energy') + return C + def _calculate_tf_ir(self, inputSignal, outputSignal, method='linear', winType=None, winSize=None, overlap=None, regularization=False): @@ -1405,41 +1456,26 @@ def _calculate_tf_ir(self, inputSignal, outputSignal, method='linear', elif inputSignal.samplingRate != outputSignal.samplingRate: raise ValueError("Both signal-like objects must have the same\ sampling rate.") + if method == 'linear': if regularization: - data = _make_pk_spectra(inputSignal.freqSignal) - outputFreqSignal = _make_pk_spectra(outputSignal.freqSignal) - freqVector = inputSignal.freqVector - b = data * 0 + 10**(-200/20) # inside signal freq range - a = data * 0 + 1 # outside signal freq range - minFreq = np.max([inputSignal.freqMin, outputSignal.freqMin]) - maxFreq = np.min([inputSignal.freqMax, outputSignal.freqMax]) - # Calculate epsilon - eps = self._crossfade_spectruns(a, b, - [minFreq/np.sqrt(2), - minFreq], - freqVector) - if maxFreq < np.min([maxFreq*np.sqrt(2), - inputSignal.samplingRate/2]): - eps = self._crossfade_spectruns(eps, a, - [maxFreq, - maxFreq*np.sqrt(2)], - freqVector) - eps = \ - eps \ - * float(np.max(np.abs(outputFreqSignal)))**2 \ - * 1/2 - C = np.conj(data) / \ - (np.conj(data)*data + eps) - C = _make_rms_spectra(C) - C = SignalObj(C, - 'freq', - inputSignal.samplingRate, - signalType='energy') + C = self._calculate_regu_spk(inputSignal, outputSignal) result = outputSignal * C else: result = outputSignal / inputSignal + elif method == 'linear_zp': # linear with zero-padding safe-guard + inputSignal_zp = inputSignal.zero_padding(num_zeros =\ + inputSignal.timeSignal.shape[0]) + outputSignal_zp = outputSignal.zero_padding(num_zeros =\ + outputSignal.timeSignal.shape[0]) + + if regularization: + C = self._calculate_regu_spk(inputSignal_zp, outputSignal_zp) + result = outputSignal_zp * C + else: + result = outputSignal_zp / inputSignal_zp + elif method == 'H1': if winType is None: winType = 'hann' @@ -1594,7 +1630,7 @@ def _crossfade_spectruns(self, a, b, freqLims, freqVector): f1idx = np.where(freqVector <= f1)[0][-1] totalSamples = a.shape[0] xsamples = f1idx - f0idx - win = ss.hanning(2*xsamples) + win = ss.hann(2*xsamples) rightWin = win[xsamples-1:-1] fullRightWin = np.concatenate((np.ones(f0idx), diff --git a/pytta/generate.py b/pytta/generate.py index b65e615..1979f54 100644 --- a/pytta/generate.py +++ b/pytta/generate.py @@ -266,8 +266,8 @@ def __do_sweep_windowing(inputSweep, # exact sample where the chirp reaches freqMax [Hz] freqMaxSample = np.where(freqSweep <= freqMax) freqMaxSample = len(freqSweep) - freqMaxSample[-1][-1] - windowStart = ss.hanning(2*freqMinSample) - windowEnd = ss.hanning(2*freqMaxSample) + windowStart = ss.hann(2*freqMinSample) + windowEnd = ss.hann(2*freqMaxSample) # Uses first half of windowStart, last half of windowEnd, and a vector of # ones with the remaining length, in between the half windows @@ -459,7 +459,7 @@ def __do_noise_windowing(inputNoise, window): # sample equivalent to the first five percent of noise duration fivePercentSample = int((5/100) * (noiseSamples)) - windowStart = ss.hanning(2*fivePercentSample) + windowStart = ss.hann(2*fivePercentSample) fullWindow = np.concatenate((windowStart[0:fivePercentSample], np.ones(int(noiseSamples-fivePercentSample)))) newNoise = (fullWindow * inputNoise.T).T