From 0a3c8dc994ebfa1b14ce4c0673b4a879571ebf8d Mon Sep 17 00:00:00 2001 From: EduardoGFilho <121968920+EduardoGFilho@users.noreply.github.com> Date: Fri, 14 Apr 2023 17:04:19 -0300 Subject: [PATCH 1/2] Fixed ak_dtfmtx The function was using three if statements, which caused the else portion to run every time the last if wasn't attended, replacing the 'if's with 'elif's fixes that. --- lasse/transforms/dft.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lasse/transforms/dft.py b/lasse/transforms/dft.py index c2b12e5..d15af43 100644 --- a/lasse/transforms/dft.py +++ b/lasse/transforms/dft.py @@ -26,10 +26,10 @@ def ak_fftmtx(N, option=1): if option == 1: # orthonormal (also called unitary) Ah = Ah / np.sqrt(N) A = np.conj(Ah) - if option == 2: # read X(k) in Volts in case x(n) is in Volts + elif option == 2: # read X(k) in Volts in case x(n) is in Volts A = np.conj(Ah) Ah = Ah / N - if option == 3: # as in Matlab/Octave: Ah = Ah + elif option == 3: # as in Matlab/Octave: Ah = Ah A = np.conj(Ah) / N else: raise Exception("Invalid option value: " + str(option)) From 4dec9d22e03663425dfb08ca42ea951f6e440a14 Mon Sep 17 00:00:00 2001 From: Eduardo Filho Date: Mon, 5 Feb 2024 08:42:22 -0300 Subject: [PATCH 2/2] rectpuls function and porting of some ak_... functions --- lasse/util/ak_functions.py | 311 +++++++++++++++++++++++++++++++++ lasse/util/python_as_matlab.py | 11 +- 2 files changed, 321 insertions(+), 1 deletion(-) create mode 100644 lasse/util/ak_functions.py diff --git a/lasse/util/ak_functions.py b/lasse/util/ak_functions.py new file mode 100644 index 0000000..9db1305 --- /dev/null +++ b/lasse/util/ak_functions.py @@ -0,0 +1,311 @@ +""" +Python reimplementation of professor Aldebaro's ak_... matlab functions +""" +# Though primary functionality has been transferred, some checks +# and customization parameters are yet to be implemented + +import numpy as np +import matplotlib.pyplot as plt +from scipy.signal import welch, get_window + + +# function bellow is adapted from wiki.python.org/moin/NumericAndScientificRecipes +""" +Find 2^n that is equal to or greater than. +""" +def nextpow2(i): + n = 0 + while 2**n < i: n +=1 + return n + +def ak_psd(x,Fs): + #Returns the PSD SdBmHz of x in dBm/Hz in the frequency range + #f=[-Fs/2,Fs/2[. Usage example: + #Fs=100 + #x=np.exp(j*2*pi*7/Fs*np.arange(40000+1)) + #S,f=ak_psd(x,Fs) + #plt.plot(f,S) #see peak at 7 Hz + N=len(x) + if N > 2048 : + M = np.round(len(x)/8).astype(int) #number of samples per block for Welch's + else: #use a single FFT in case N is not longer than 2048 + M=len(x) + + b=nextpow2(M) + Nfft=2**b #choose a power of 2 for FFT-length + ham_window = get_window('hamming',M, fftbins = False) + _,S = welch(x = x, fs = Fs, window = ham_window,nfft=Nfft, return_onesided = False) + df=Fs/Nfft #FFT frequency spacing + f=np.arange(-Fs/2,Fs/2,df) - df #frequency (abscissa) axis + S=np.fft.fftshift(S) #move negative part to the left + S=S*1000 #convert from Watts to mWatts + SdBmHz = 10*np.log10(S) #provide PSD in dBm/Hz + #plt.plot(f,SdBmHz) + return SdBmHz,f + +def ak_sinc_reconstruction(n, x, Ts, n_oversampled, xo, textra = 0): + + #function [x_reconstructed, x_parcels, t_oversampled, ... + # t_oversampled_expanded] = ak_sinc_reconstruction(n, x, Ts, ... + # n_oversampled, xo, textra) + ## Illustrate perfect sinc reconstruction + #Inputs: + #n - vector with integers as abscissa of x + #x - samples of signal to be reconstructed + #Ts - sampling interval (in seconds) + #n_oversampled - vector with integers as abscissa of oversampled x + #xo - oversampled x, to be used for comparison purposes + #textra - extra time (in s) to be added for visualization purposes + #Outputs: + #x_reconstructed - reconstructed signal + #x_parcels - matrix with individual sincs used in x_reconstructed + #t_oversampled - discrete-time in seconds, corresponding to + # input n_oversampled multiplied by the oversampled Ts + #t_oversampled_expanded - similar to t_oversampled but with the + # # addition of textra in the beginning and end + + # TODO: Implement checks in python + # if nargin < 6 + # textra = 0 #5*Ts + # end + # if sum(logical(rem(n,1))) ~= 0 #check if vector has only integers + # error('Elements of vector n must be integers!') + # end + # if sum(logical(rem(n_oversampled,1))) ~= 0 #check if integers + # error('Elements of vector n_oversampled must be integers!') + # end + + num_discrete_samples = len(x) + # if length(n) ~= num_discrete_samples + # error('Lengths of n and xn must be the same!') + # end + + num_oversampled_samples = len(xo) + # if length(n_oversampled) ~= num_oversampled_samples + # error('Lengths of n_oversampled and xo must be the same!') + # end + + #oversampling fac is adapted from wiki.python.org/mor can be found by + #N_o = L (N-1) + 1, where N is the number of samples + #in x and L is the oversampling factor + oversampling_factor = (num_oversampled_samples-1) / (num_discrete_samples-1) + # if mod(oversampling_factor,1) ~= 0 + # error('Oversampling factor must be an integer!') + # end + oversampled_Ts = Ts/oversampling_factor + + ## Create expanded time axis to help visualation if requested + #get number of samples to be added before and after original signal + nextra = np.ceil(textra/oversampled_Ts).astype(int) + + n_oversampled_expanded = np.arange(n_oversampled[0] - nextra, n_oversampled[-1] + nextra) + t_oversampled_expanded = n_oversampled_expanded*oversampled_Ts + + total_samples = len(t_oversampled_expanded) + ## Find parcel corresponding to each sinc and sum them to + # compose x_reconstructed + x_parcels = np.zeros((num_discrete_samples, total_samples)) + x_reconstructed = np.zeros(total_samples) + for ncounter in range(num_discrete_samples): + this_n = n[ncounter] #value of n, for instance, n=-4 + #Sinc delayed by this_n*Ts: + this_sinc=x[ncounter]*np.sinc((t_oversampled_expanded-this_n*Ts)/Ts) + x_parcels[ncounter,:]=this_sinc #save this sinc + x_reconstructed = x_reconstructed + this_sinc #add sinc contribution + + #create vector in case one wants to plot results + t_oversampled = n_oversampled*oversampled_Ts + # if nargout < 1 + # #plot results + # ak_plot_sinc_reconstruction(n, x, Ts, t_oversampled, xo, ... + # t_oversampled_expanded,x_reconstructed,x_parcels) + # end + + ak_plot_sinc_reconstruction(n, x, Ts, t_oversampled, xo,\ + t_oversampled_expanded,x_reconstructed,x_parcels) + + return x_reconstructed, x_parcels, t_oversampled, t_oversampled_expanded + +def ak_plot_sinc_reconstruction(n, x, Ts, t_oversampled, xo,\ + t_oversampled_expanded,x_reconstructed,x_parcels): + #function ak_plot_sinc_reconstruction(n, x, Ts, t_oversampled, xo,... + # t_oversampled_expanded,x_reconstructed,x_parcels) + #See ak_sinc_reconstruction.m + + plt.figure(1) + plt.subplot(411) + plt.plot(t_oversampled, xo) + plt.xlabel('t (s)') + plt.ylabel('Original x(t)') + plt.tight_layout() + + plt.subplot(413) + plt.stem(n, x) + plt.xlabel('Discrete-time n') + plt.ylabel('x[n]') + plt.grid() + plt.tight_layout() + + plt.subplot(412) + t = n*Ts + #T=0.2; ak_sampledsignalsplot(x,[],T,'color','r') + ak_sampledsignalsplot(x, t, []) + + #plt.plot(t, np.zeros(len(t)),color = 'b') + plt.xlabel('t (s)') + plt.ylabel('x$_s$(t)') + plt.tight_layout() + plt.grid() + + plt.subplot(414) + plt.plot(t_oversampled, xo) + + plt.plot(t_oversampled_expanded,x_reconstructed,'r-'); #,'--','LineWidth',1.5); + + xmin, xmax, ymin, ymax = plt.axis() + support = t_oversampled[-1] - t_oversampled[0] + xmin = t_oversampled[0] - support*0.05 + xmax = t_oversampled[-1] + support*0.05 + + plt.axis([xmin, xmax, ymin, ymax]) + plt.legend(['Original','Reconstructed']) + plt.tight_layout() + plt.ylabel('x(t)') + plt.xlabel('t (s)') + + #ak_changeFigureSize(1.5, 1.5) #expand figure + + ## Show the parcels (each sinc) in sinc reconstruction + plt.figure(2) + + plt.plot(t_oversampled_expanded,x_reconstructed, linewidth = 1.5) + + # TODO: implement colors + #C=colororder; #Octave does not recognize it + # C=get(gca,"colororder"); + # newcolors = [0.83 0.14 0.14 + # 1.00 0.54 0.00 + # 0.47 0.25 0.80]; + # C=[C(1,:) + # newcolors + # C(2:end,:)]; + # #colororder(C); #Octave does not recognize it + # set(gca,"colororder",C); + # [num_colors, ~] = size(C); + # hold on + + num_discrete_samples = len(x) + for i in range(num_discrete_samples): + plt.plot(t_oversampled_expanded,x_parcels[i,:],'--') + plt.stem(n*Ts,x) + plt.plot(n[i]*Ts,x[i], marker = 'x', markersize = 20, linewidth = 2.5)#,'x','MarkerSize',20,'Color',C(j,:),'LineWidth',2.5); + + y_zeros = np.arange(t_oversampled_expanded[0],t_oversampled_expanded[-1]+1,Ts) + plt.stem(y_zeros , np.zeros(len(y_zeros)), markerfmt = 'ko') + plt.xlabel('t (s)') + plt.ylabel('Reconstructed x(t)') + + plt.tight_layout() + #ak_changeFigureSize(1.5, 1.5) #expand figure + +def ak_changeFigureSize(width_factor,height_factor): + # function ak_changeFigureSize(width_factor, height_factor) + #Change figure size by the factors width_factor and height_factor, + #and position figure in center of screen. + + #From https://www.cs.cornell.edu/courses/cs100m/2007fa/Graphics/Position.pdf + #set(gcf,'position',[a b H W]) + #places the lower left corner of an H-by-W figure window at (a, b). + + w,h = plt.gcf().get_size_inches() + plt.gcf().set_size_inches(w*width_factor,h*height_factor) + # pos=get(gcf, 'Position'); #get figure's position on screen + # pos(3)=floor(pos(3)*width_factor); #adjust the width + # pos(4)=floor(pos(4)*height_factor); #adjust the height + # set(gcf,'Position',pos); + # movegui(gcf,'center') + +# TODO: implement 'varargin' +def ak_sampledsignalsplot(x,t,T): + # function ak_sampledsignalsplot(x,t,T,varargin) + #Plots the time series x as impulses. The abscissa is t as time-axis or + #generated automatically with T as sampling interval. + #The distinction with respect to ak_impulseplot is that this function does + #not plot the zero values + #Examples: + # x=1:3; t=1:0.5:2; ak_sampledsignalsplot(x,t,[],'color','r') + # T=0.5; ak_sampledsignalsplot(x,[],T,'color','r') + + # TODO: implement syntax check + #Check syntax: + # if isempty(t) + # if isempty(T) + # error(['You need to specify the time-axis t or the ' ... + # 'sampling interval T']); + # end + # if ~isnumeric(T) + # error('Syntax error: T must be numeric (when specified)'); + # end + # #generate time-axis + # t=0:T:(length(x)-1)*T; + # else + # if nargin > 2 #allow the user specify only two parameters + # if ~isempty(T) + # error(['You cannot specify both the time-axis t and the ' ... + # 'sampling interval T. Check the syntax']); + # end + # end + # end + + #AK April 24 - 2021 + #because negative impulses were having different colors than positive ones + #I am defining a default color: + # if isempty(varargin) + # varargin={'color','b'}; + # end + + #remember whether hold is on or off + # holdison = 0; + # if ishold + # holdison = 1; + # end + #add some extra comments to stem for plotting an arrowhead + #do not forget to pass it the other parameters in varargin{:} + + indicesPos = np.where(x>0) + indicesNeg = np.where(x<0) + indicesZero = np.where(x==0) + + markersize = 9 + linewidth = 2 + + pos_markerline, pos_stemline, pos_baseline = plt.stem(t[indicesPos],x[indicesPos]) + plt.setp(pos_markerline, marker = '^', markersize = markersize,\ + color = 'b') + plt.setp(pos_stemline, color = 'b', linewidth = linewidth) + plt.setp(pos_baseline,linestyle = '') + + neg_markerline, neg_stemline, neg_baseline = plt.stem(t[indicesNeg],x[indicesNeg]) + plt.setp(neg_markerline, marker = 'v', markersize = markersize,\ + color = 'b') + plt.setp(neg_stemline, color = 'b', linewidth = linewidth) + plt.setp(neg_baseline,linestyle = '') + + # if 0 #do not plot + # stem(t(indicesZero),x(indicesZero),'marker','o','markersize',markersize,... + # 'markerfacecolor','auto', ... + # 'LineStyle','-', 'LineWidth',linewidth, ... + # varargin{:}); + # end + + plt.plot(t,np.zeros(np.size(t)), color = 'b', linestyle='-') #plot line to represent zeros + + + #restore previous hold situation + # hold off + # if holdison==1 + # hold on + # end + + plt.xlabel('Time (s)') + plt.ylabel('Amplitude') \ No newline at end of file diff --git a/lasse/util/python_as_matlab.py b/lasse/util/python_as_matlab.py index c3a4a5e..e7428ef 100644 --- a/lasse/util/python_as_matlab.py +++ b/lasse/util/python_as_matlab.py @@ -94,6 +94,15 @@ def test_xcorr(): print(Rxy) print(lags) +def rectpuls(x,w = 1): + # TODO: use warnings rather then prints + if not isinstance(x, np.ndarray): + print("x must be a NumPy array.") + # pulse with duration [-0.5,0.5) + y = 1.0*(x>0.5*w) - 1.0*(x<=-0.5*w) + y = np.array(y,dtype = x.dtype) + return y + if __name__ == "__main__": - test_xcorr() + test_xcorr() \ No newline at end of file