[docs]deffft1d(x:np.ndarray,y:np.ndarray,shift:bool=True)->tuple[np.ndarray,np.ndarray]:"""Compute FFT on X,Y data. Args: x: X data y: Y data shift: Shift the zero frequency to the center of the spectrum. Defaults to True. Returns: X data, Y data (tuple) """y1=np.fft.fft(y)x1=np.fft.fftfreq(x.size,d=x[1]-x[0])ifshift:x1=np.fft.fftshift(x1)y1=np.fft.fftshift(y1)returnx1,y1
[docs]defifft1d(x:np.ndarray,y:np.ndarray,shift:bool=True)->tuple[np.ndarray,np.ndarray]:"""Compute iFFT on X,Y data. Args: x: X data y: Y data shift: Shift the zero frequency to the center of the spectrum. Defaults to True. Returns: X data, Y data (tuple) """ifshift:y=np.fft.ifftshift(y)y1=np.fft.ifft(y)# Recalculate the original time domain arraydt=1.0/(x[-1]-x[0]+(x[1]-x[0]))x1=np.arange(y1.size)*dtreturnx1,y1.real
[docs]defmagnitude_spectrum(x:np.ndarray,y:np.ndarray,log_scale:bool=False)->tuple[np.ndarray,np.ndarray]:"""Compute magnitude spectrum. Args: x: X data y: Y data log_scale: Use log scale. Defaults to False. Returns: Magnitude spectrum (X data, Y data) """x1,y1=fft1d(x,y)iflog_scale:y_mag=20*np.log10(np.abs(y1))y_mag=np.abs(y1)returnx1,y_mag
[docs]defphase_spectrum(x:np.ndarray,y:np.ndarray)->tuple[np.ndarray,np.ndarray]:"""Compute phase spectrum. Args: x: X data y: Y data Returns: Phase spectrum in degrees (X data, Y data) """x1,y1=fft1d(x,y)y_phase=np.rad2deg(np.angle(y1))returnx1,y_phase
[docs]defpsd(x:np.ndarray,y:np.ndarray,log_scale:bool=False)->tuple[np.ndarray,np.ndarray]:"""Compute Power Spectral Density (PSD), using the Welch method. Args: x: X data y: Y data log_scale: Use log scale. Defaults to False. Returns: Power Spectral Density (PSD): X data, Y data (tuple) """x1,y1=scipy.signal.welch(y,fs=sampling_rate(x))iflog_scale:y1=10*np.log10(y1)returnx1,y1
[docs]defsort_frequencies(x:np.ndarray,y:np.ndarray)->np.ndarray:"""Sort from X,Y data by computing FFT(y). Args: x: X data y: Y data Returns: Sorted frequencies in ascending order """freqs,fourier=fft1d(x,y,shift=False)returnfreqs[np.argsort(fourier)]
[docs]defpeak_indices(y,thres:float=0.3,min_dist:int=1,thres_abs:bool=False)->np.ndarray:# Copyright (c) 2014 Lucas Hermann Negri# Unmodified code snippet from PeakUtils 1.3.0"""Peak detection routine. Finds the numeric index of the peaks in *y* by taking its first order difference. By using *thres* and *min_dist* parameters, it is possible to reduce the number of detected peaks. *y* must be signed. Parameters ---------- y : ndarray (signed) 1D amplitude data to search for peaks. thres : float between [0., 1.] Normalized threshold. Only the peaks with amplitude higher than the threshold will be detected. min_dist : int Minimum distance between each detected peak. The peak with the highest amplitude is preferred to satisfy this constraint. thres_abs: boolean If True, the thres value will be interpreted as an absolute value, instead of a normalized threshold. Returns ------- ndarray Array containing the numeric indices of the peaks that were detected """ifisinstance(y,np.ndarray)andnp.issubdtype(y.dtype,np.unsignedinteger):raiseValueError("y must be signed")ifnotthres_abs:thres=thres*(np.max(y)-np.min(y))+np.min(y)# compute first order differencedy=np.diff(y)# propagate left and right values successively to fill all plateau pixels# (0-value)(zeros,)=np.where(dy==0)# check if the signal is totally flatiflen(zeros)==len(y)-1:returnnp.array([])iflen(zeros):# compute first order difference of zero indiceszeros_diff=np.diff(zeros)# check when zeros are not chained together(zeros_diff_not_one,)=np.add(np.where(zeros_diff!=1),1)# make an array of the chained zero indiceszero_plateaus=np.split(zeros,zeros_diff_not_one)# fix if leftmost value in dy is zeroifzero_plateaus[0][0]==0:dy[zero_plateaus[0]]=dy[zero_plateaus[0][-1]+1]zero_plateaus.pop(0)# fix if rightmost value of dy is zeroiflen(zero_plateaus)>0andzero_plateaus[-1][-1]==len(dy)-1:dy[zero_plateaus[-1]]=dy[zero_plateaus[-1][0]-1]zero_plateaus.pop(-1)# for each chain of zero indicesforplateauinzero_plateaus:median=np.median(plateau)# set leftmost values to leftmost non zero valuesdy[plateau[plateau<median]]=dy[plateau[0]-1]# set rightmost and middle values to rightmost non zero valuesdy[plateau[plateau>=median]]=dy[plateau[-1]+1]# find the peaks by using the first order differencepeaks=np.where((np.hstack([dy,0.0])<0.0)&(np.hstack([0.0,dy])>0.0)&(np.greater(y,thres)))[0]# handle multiple peaks, respecting the minimum distanceifpeaks.size>1andmin_dist>1:highest=peaks[np.argsort(y[peaks])][::-1]rem=np.ones(y.size,dtype=bool)rem[peaks]=Falseforpeakinhighest:ifnotrem[peak]:sl=slice(max(0,peak-min_dist),peak+min_dist+1)rem[sl]=Truerem[peak]=Falsepeaks=np.arange(y.size)[~rem]returnpeaks
[docs]defxpeak(x:np.ndarray,y:np.ndarray)->float:"""Return default peak X-position (assuming a single peak). Args: x: X data y: Y data Returns: Peak X-position """peaks=peak_indices(y)ifpeaks.size==1:returnx[peaks[0]]returnnp.average(x,weights=y)
[docs]definterpolate(x:np.ndarray,y:np.ndarray,xnew:np.ndarray,method:Literal["linear","spline","quadratic","cubic","barycentric","pchip"],fill_value:float|None=None,)->np.ndarray:"""Interpolate data. Args: x: X data y: Y data xnew: New X data method: Interpolation method fill_value: Fill value. Defaults to None. This value is used to fill in for requested points outside of the X data range. It is only used if the method argument is 'linear', 'cubic' or 'pchip'. Returns: Interpolated Y data """interpolator_extrap=Noneifmethod=="linear":# Linear interpolation using NumPy's interp function:ynew=np.interp(xnew,x,y,left=fill_value,right=fill_value)elifmethod=="spline":# Spline using 1-D interpolation with SciPy's interpolate package:# pylint: disable=unbalanced-tuple-unpackingknots,coeffs,degree=scipy.interpolate.splrep(x,y,s=0)ynew=scipy.interpolate.splev(xnew,(knots,coeffs,degree),der=0)elifmethod=="quadratic":# Quadratic interpolation using NumPy's polyval function:coeffs=np.polyfit(x,y,2)ynew=np.polyval(coeffs,xnew)elifmethod=="cubic":# Cubic interpolation using SciPy's Akima1DInterpolator class:interpolator_extrap=scipy.interpolate.Akima1DInterpolator(x,y)elifmethod=="barycentric":# Barycentric interpolation using SciPy's BarycentricInterpolator class:interpolator=scipy.interpolate.BarycentricInterpolator(x,y)ynew=interpolator(xnew)elifmethod=="pchip":# PCHIP interpolation using SciPy's PchipInterpolator class:interpolator_extrap=scipy.interpolate.PchipInterpolator(x,y)else:raiseValueError(f"Invalid interpolation method {method}")ifinterpolator_extrapisnotNone:ynew=interpolator_extrap(xnew,extrapolate=fill_valueisNone)iffill_valueisnotNone:ynew[xnew<x[0]]=fill_valueynew[xnew>x[-1]]=fill_valuereturnynew
[docs]defwindowing(y:np.ndarray,method:Literal["barthann","bartlett","blackman","blackman-harris","bohman","boxcar","cosine","exponential","flat-top","hamming","hanning","lanczos","nuttall","parzen","rectangular","taylor","tukey","kaiser","gaussian",]="hamming",alpha:float=0.5,beta:float=14.0,sigma:float=7.0,)->np.ndarray:"""Apply windowing to the input data. Args: x: X data y: Y data method: Windowing function. Defaults to "hamming". alpha: Tukey window parameter. Defaults to 0.5. beta: Kaiser window parameter. Defaults to 14.0. sigma: Gaussian window parameter. Defaults to 7.0. Returns: Windowed Y data """# Cases without parameters:win_func={"barthann":scipy.signal.windows.barthann,"bartlett":np.bartlett,"blackman":np.blackman,"blackman-harris":scipy.signal.windows.blackmanharris,"bohman":scipy.signal.windows.bohman,"boxcar":scipy.signal.windows.boxcar,"cosine":scipy.signal.windows.cosine,"exponential":scipy.signal.windows.exponential,"flat-top":scipy.signal.windows.flattop,"hamming":np.hamming,"hanning":np.hanning,"lanczos":scipy.signal.windows.lanczos,"nuttall":scipy.signal.windows.nuttall,"parzen":scipy.signal.windows.parzen,"rectangular":np.ones,"taylor":scipy.signal.windows.taylor,}.get(method)ifwin_funcisnotNone:returny*win_func(len(y))# Cases with parameters:ifmethod=="tukey":returny*scipy.signal.windows.tukey(len(y),alpha)ifmethod=="kaiser":returny*np.kaiser(len(y),beta)ifmethod=="gaussian":returny*scipy.signal.windows.gaussian(len(y),sigma)raiseValueError(f"Invalid window type {method}")
[docs]@classmethoddefget_amp_from_amplitude(cls,amplitude,sigma):"""Return amp from function amplitude and sigma"""returnamplitude*(sigma*np.sqrt(2*np.pi))
[docs]@classmethoddefamplitude(cls,amp,sigma):"""Return function amplitude"""returnamp/(sigma*np.sqrt(2*np.pi))
[docs]@classmethoddeffwhm(cls,amp,sigma):"""Return function FWHM"""return2*sigma*np.sqrt(2*np.log(2))
[docs]classLorentzianModel(FitModel):"""1-dimensional Lorentzian fit model"""
[docs]@classmethoddeffwhm(cls,amp,sigma):"""Return function FWHM"""wg=GaussianModel.fwhm(amp,sigma)wl=LorentzianModel.fwhm(amp,sigma)return0.5346*wl+np.sqrt(0.2166*wl**2+wg**2)
[docs]deffind_nearest_zero_point_idx(y:np.ndarray)->np.ndarray:"""Find the x indices where the corresponding y is the closest to zero Args: y: Y data Returns: Indices of the points right before or at zero crossing """xi=np.where((y[:-1]>=0)&(y[1:]<=0)|(y[:-1]<=0)&(y[1:]>=0))[0]returnxi
[docs]deffind_x_at_value(x:np.ndarray,y:np.ndarray,value:float)->np.ndarray:"""Find the x value where the y value is the closest to the given value using linear interpolation to deduce the precise x value. Args: x: X data y: Y data value: Value to find Returns: X value where the Y value is the closest to the given value """leveled_y=y-valuexi_before=find_nearest_zero_point_idx(leveled_y)xi_after=xi_before+1iflen(xi_before)==0:returnnp.array([0.0])# linear interpolationp=(leveled_y[xi_after]-leveled_y[xi_before])/(x[xi_after]-x[xi_before])ori=leveled_y[xi_after]-p*x[xi_after]x0=-ori/p# where the curve cut the absissareturnx0
[docs]defbandwidth(data:np.ndarray,level:float=3.0)->float:"""Compute the bandwidth of the signal at a given level. Args: data: X,Y data level: Level in dB at which the bandwidth is computed. Defaults to 3.0. Returns: Bandwidth of the signal at the given level: segment coordinates """x,y=datahalf_max:float=np.max(y)-levelbw=find_x_at_value(x,y,half_max)[0]coords=(x[0],half_max,bw,half_max)returncoords
[docs]defsinusoidal_model(x:np.ndarray,a:float,f:float,phi:float,offset:float)->np.ndarray:"""Sinusoidal model function."""returna*np.sin(2*np.pi*f*x+phi)+offset
[docs]defsinusoidal_fit(x:np.ndarray,y:np.ndarray)->tuple[tuple[float,float,float,float],float]:"""Fit a sinusoidal model to the input data. Args: x: X data y: Y data Returns: A tuple containing the fit parameters (amplitude, frequency, phase, offset) and the residuals """# Initial guess for the parameters# ==================================================================================offset=np.mean(y)amp=(np.max(y)-np.min(y))/2phase_origin=0# Search for the maximum of the FFTi_maxfft=np.argmax(np.abs(np.fft.fft(y-offset)))ifi_maxfft>len(x)/2:# If the index is greater than N/2, we are in the mirrored half spectrum# (negative frequencies)i_maxfft=len(x)-i_maxfftfreq=i_maxfft/(x[-1]-x[0])# ==================================================================================defoptfunc(fitparams:np.ndarray,x:np.ndarray,y:np.ndarray)->np.ndarray:"""Optimization function."""returny-sinusoidal_model(x,*fitparams)# Fit the model to the datafitparams=scipy.optimize.leastsq(optfunc,[amp,freq,phase_origin,offset],args=(x,y))[0]y_th=sinusoidal_model(x,*fitparams)residuals=np.std(y-y_th)returnfitparams,residuals
[docs]defsinus_frequency(x:np.ndarray,y:np.ndarray)->float:"""Compute the frequency of a sinusoidal signal. Args: x: x signal data y: y signal data Returns: Frequency of the sinusoidal signal """fitparams,_residuals=sinusoidal_fit(x,y)returnfitparams[1]
[docs]defenob(x:np.ndarray,y:np.ndarray,full_scale:float=1.0)->float:"""Compute Effective Number of Bits (ENOB). Args: x: x signal data y: y signal data full_scale: Full scale(V). Defaults to 1.0. Returns: Effective Number of Bits (ENOB) """_fitparams,residuals=sinusoidal_fit(x,y)return-np.log2(residuals*np.sqrt(12)/full_scale)
[docs]defsinad(x:np.ndarray,y:np.ndarray,full_scale:float=1.0,unit:Literal["dBc","dBFS"]="dBc",)->float:"""Compute Signal-to-Noise and Distortion Ratio (SINAD). Args: x: x signal data y: y signal data full_scale: Full scale(V). Defaults to 1.0. unit: Unit of the input data. Valid values are 'dBc' and 'dBFS'. Defaults to 'dBc'. Returns: Signal-to-Noise and Distortion Ratio (SINAD) """fitparams,residuals=sinusoidal_fit(x,y)amp=fitparams[0]# Compute the power of the fundamentalpowf=np.abs(amp/np.sqrt(2))ifunit=="dBc"elsefull_scale/(2*np.sqrt(2))return20*np.log10(powf/residuals)
[docs]defthd(x:np.ndarray,y:np.ndarray,full_scale:float=1.0,unit:Literal["dBc","dBFS"]="dBc",nb_harm:int=5,)->float:"""Compute Total Harmonic Distortion (THD). Args: x: x signal data y: y signal data full_scale: Full scale(V). Defaults to 1.0. unit: Unit of the input data. Valid values are 'dBc' and 'dBFS'. Defaults to 'dBc'. nb_harm: Number of harmonics to consider. Defaults to 5. Returns: Total Harmonic Distortion (THD) """fitparams,_residuals=sinusoidal_fit(x,y)offset=np.mean(y)amp,freq=fitparams[:2]ampfft=np.abs(np.fft.fft(y-offset))# Compute the power of the fundamentalifunit=="dBc":powfund=np.max(ampfft[:len(ampfft)//2])else:powfund=(full_scale/(2*np.sqrt(2)))*(len(x)/np.sqrt(2))sumharm=0foriinnp.arange(nb_harm+2)[2:]:a=i*np.ceil(freq*(x[-1]-x[0]))amp=ampfft[int(a-5):int(a+5)]iflen(amp)>0:sumharm+=np.max(amp)return20*np.log10(sumharm/powfund)
[docs]defsfdr(x:np.ndarray,y:np.ndarray,full_scale:float=1.0,unit:Literal["dBc","dBFS"]="dBc",)->float:"""Compute Spurious-Free Dynamic Range (SFDR). Args: x: x signal data y: y signal data full_scale: Full scale(V). Defaults to 1.0. unit: Unit of the input data. Valid values are 'dBc' and 'dBFS'. Defaults to 'dBc'. Returns: Spurious-Free Dynamic Range (SFDR) """fitparams,_residuals=sinusoidal_fit(x,y)# Compute the power of the fundamentalifunit=="dBc":powfund=np.max(np.abs(np.fft.fft(y)))else:powfund=(full_scale/(2*np.sqrt(2)))*(len(x)/np.sqrt(2))maxspike=np.max(np.abs(np.fft.fft(y-sinusoidal_model(x,*fitparams))))return20*np.log10(powfund/maxspike)
[docs]defsnr(x:np.ndarray,y:np.ndarray,full_scale:float=1.0,unit:Literal["dBc","dBFS"]="dBc",)->float:"""Compute Signal-to-Noise Ratio (SNR). Args: x: x signal data y: y signal data full_scale: Full scale(V). Defaults to 1.0. unit: Unit of the input data. Valid values are 'dBc' and 'dBFS'. Defaults to 'dBc'. Returns: Signal-to-Noise Ratio (SNR) """fitparams,_residuals=sinusoidal_fit(x,y)# Compute the power of the fundamentalifunit=="dBc":powfund=np.max(np.abs(np.fft.fft(y)))else:powfund=(full_scale/(2*np.sqrt(2)))*(len(x)/np.sqrt(2))noise=np.sqrt(np.mean((y-sinusoidal_model(x,*fitparams))**2))return20*np.log10(powfund/noise)
[docs]defsampling_period(x:np.ndarray)->float:"""Compute sampling period Args: x: X data Returns: Sampling period """steps=np.diff(x)ifnotnp.isclose(np.diff(steps).max(),0,atol=1e-10):warnings.warn("Non-constant sampling signal")returnsteps[0]
[docs]defsampling_rate(x:np.ndarray)->float:"""Compute mean sampling rate Args: x: X data Returns: Sampling rate """return1.0/sampling_period(x)
[docs]deffwhm(data:np.ndarray,method:Literal["zero-crossing","gauss","lorentz","voigt"]="zero-crossing",xmin:float|None=None,xmax:float|None=None,)->tuple[float,float,float,float]:"""Compute Full Width at Half Maximum (FWHM) of the input data Args: data: X,Y data method: Calculation method. Two types of methods are supported: a zero-crossing method and fitting methods (based on various models: Gauss, Lorentz, Voigt). Defaults to "zero-crossing". xmin: Lower X bound for the fitting. Defaults to None (no lower bound, i.e. the fitting starts from the first point). xmax: Upper X bound for the fitting. Defaults to None (no upper bound, i.e. the fitting ends at the last point) Returns: FWHM segment coordinates """x,y=datadx,dy,base=np.max(x)-np.min(x),np.max(y)-np.min(y),np.min(y)sigma,mu=dx*0.1,xpeak(x,y)ifisinstance(xmin,float):indices=np.where(x>=xmin)[0]x=x[indices]y=y[indices]ifisinstance(xmax,float):indices=np.where(x<=xmax)[0]x=x[indices]y=y[indices]ifmethod=="zero-crossing":hmax=dy*0.5+np.min(y)fx=find_x_at_value(x,y,hmax)iffx.size!=2:warnings.warn(f"Ambiguous zero-crossing points (found {fx.size} points)")returnfx[0],hmax,fx[-1],hmaxtry:FitModelClass:type[FitModel]={"gauss":GaussianModel,"lorentz":LorentzianModel,"voigt":VoigtModel,}[method]exceptKeyErrorasexc:raiseValueError(f"Invalid method {method}")fromexcdeffunc(params):"""Fitting model function"""# pylint: disable=cell-var-from-loopreturny-FitModelClass.func(x,*params)amp=FitModelClass.get_amp_from_amplitude(dy,sigma)(amp,sigma,mu,base),_ier=scipy.optimize.leastsq(func,np.array([amp,sigma,mu,base]))returnFitModelClass.half_max_segment(amp,sigma,mu,base)
[docs]deffw1e2(data:np.ndarray)->tuple[float,float,float,float]:"""Compute Full Width at 1/e� of the input data (using a Gaussian model fitting). Args: data: X,Y data Returns: FW at 1/e² segment coordinates """x,y=datadx,dy,base=np.max(x)-np.min(x),np.max(y)-np.min(y),np.min(y)sigma,mu=dx*0.1,xpeak(x,y)amp=GaussianModel.get_amp_from_amplitude(dy,sigma)p_in=np.array([amp,sigma,mu,base])deffunc(params):"""Fitting model function"""# pylint: disable=cell-var-from-loopreturny-GaussianModel.func(x,*params)p_out,_ier=scipy.optimize.leastsq(func,p_in)amp,sigma,mu,base=p_outhw=2*sigmayhm=GaussianModel.amplitude(amp,sigma)/np.e**2+basereturnmu-hw,yhm,mu+hw,yhm