# -*- coding: utf-8 -*- """ Routines for signal processing. Coded by József Bór """ import matplotlib.pyplot as plt import matplotlib.dates as mdates import numpy as np import datetime as dt def Spectrogram(x, fs=1.0, nperwindow=None, windowfunc=None, nfft=None, \ slidewindowbyn=None, sgram_type=None): """ Calculates the spectrogram of the real valued time series in x Parameters: x : numpy array or list of real values fs : sampling frequency (optional, defaults to 1.0) nperwindow : number of data points for single spectrum calculation Defaults to the length of array x windowfunc : name of window function to use on each data segment numpy.window is used to calculate the window function values available windows: 'square' = no window (default) 'hanning' = Hanning window 'hamming' = Hamming window 'blackman' = Blackman window nfft : number of points the FFT is calculated from nfft = nperwindow by default if nfft > nperwindow, data in the window will be zero padded if nfft < nperwindow, data in the window will be truncated Note that window function will be applied first, extension and truncation to nfft comes second slidewindowbyn : data window will be slid by this number of points to calculate the next FFT. The FFT will be calculated only if the next data window is filled with data completely. Defaults to the length of array x sgram_type : set the output of the program 'ampl' = return err, x_indxs, f, ampl (the amplitude spectrogram) this is the default setting 'phas' = return err, x_indxs, f, phas (the phase spectrogram) 'both' = return err, x_indxs, f, ampl, phas The code returns (depending on the value of 'sgram_type' argument) err : string that contains the description of any errors or warnings Upon an error, empty lists are returned as follows err, [], [], [], [] x_indxs : array of indexes in x to which calculated spectra are assigned. The index of first data for each data window is put in this array f : frequencies ampl : 2 dimensional nd array of spectral amplitudes, if requested In ampl[i][j], i corresponds to the frequency f[i] j corresponds to the spectrum at x_indxs[j] amplitudes are calculted directly from complex coefficients returned by the numpy.rfft() function amplitudes are normalized by nperwindow. No other normalization is done and the reduction due to a window function (if applied) is not corrected phas : 2 dimensional nd array of phase spectrum (radians), if requested Phase values are in the range [-pi..pi] In phas[i][j], i corresponds to the frequency f[i] j corresponds to the spectrum x_indxs[j] See and run the function Spectrogram_test_routine() in this module for an example. The code was developed under WinPython-32bit-2.7.6.4 This function requires numpy. The test function in the modul uses also the modules datetime and matplotlib. """ #### Parse and check function arguments # sgarm_type if sgram_type == None : sgram_type = 'ampl' if sgram_type == 'ampl' : pass elif sgram_type == 'phas' : pass elif sgram_type == 'both' : pass else: err = 'Error: invalid setting for argument "sgram_type"' return err, [], [], [], [] # windowfunc if windowfunc == None : sgram_type = 'square' if windowfunc == 'square' : pass elif windowfunc == 'hanning' : pass elif windowfunc == 'hamming' : pass elif windowfunc == 'blackman' : pass else: err = 'Error: invalid setting for argument "windowfunc"' return err, [], [], [], [] #nperwindow if nperwindow == None : nperwindow = x.size if nperwindow <= 0 or nperwindow > x.size: err = 'Error: "nperwindow" must be a positive integer up to the' \ +'number of elements in "x"' return err, [], [], [], [] #nfft if nfft == None : nfft = nperwindow #slidewindowbyn if slidewindowbyn == None : slidewindowbyn = x.size if slidewindowbyn <= 0: err = 'Error: "slidewindowbyn" must be a positive integer' return err, [], [], [], [] #### Initiate error string err = '' #### Prepare the window function codestr= 'np.{0}({1})'.format(windowfunc,nperwindow) wnd = eval(codestr) #### prepare the frequency array f = np.fft.rfftfreq(nfft, 1.0/fs) #### prepare the array of x indexes x_indxs = np.arange(0, x.size, slidewindowbyn) x_indxs = x_indxs[np.where(x_indxs+nperwindow-1 < x.size)[0]] #### prepare the 2 dimensional array for fft rfft_x = np.ndarray(shape=(f.size,x_indxs.size),dtype=complex) #### calculate the complex dynamic spectrum using rfft for k in range(x_indxs.size): # extract the next segment segment = np.copy(x[x_indxs[k]:x_indxs[k]+nperwindow]) # Apply window function segment *= wnd rfft_x[:,k] = np.fft.rfft(segment,nfft) #### create the requested spectrograms if sgram_type == 'ampl' or sgram_type == 'both': ampl = np.absolute(rfft_x)/float(nperwindow) if sgram_type == 'phas' or sgram_type == 'both': phas = np.angle(rfft_x) if sgram_type == 'ampl' : return err, x_indxs, f, ampl elif sgram_type == 'phas' : return err, x_indxs, f, phas return err, x_indxs, f, ampl, phas def Chirp_T(freq_start=1.0, freq_end=10.0, duration_sec=1.0, sampling_freq=100.0, initial_phase_rad=0.0): """ Creates a linear chirp waveform and returns it as a list """ # just for testing #freq_start = 10 #freq_end = 20 #duration_sec = 2 #sampling_freq = 100 #initial_phase_rad = np.pi/3 # the time axis t = np.arange(0,duration_sec,1.0/sampling_freq) # the rate of frequency change k = (freq_end-freq_start)/duration_sec # the chirp waveform x = [np.cos(2*np.pi*(freq_start+i*k/2.0)*i + initial_phase_rad) for i in t] # just for testing #plt.plot(t,x) return x def Chirp_N(freq_start=1.0, freq_end=1.0, sampling_freq=100.0, n_pts=100, initial_phase_rad=0.0): """ Creates a linear chirp waveform and returns it as a list """ # just for testing #freq_start = 10 #freq_end = 20 #sampling_freq = 100 #n_pts = 200 #initial_phase_rad = np.pi/3 # the time axis t = [i*1.0/sampling_freq for i in np.arange(n_pts)] # the rate of frequency change k = (freq_end-freq_start)/t[-1] # the chirp waveform x = [np.cos(2*np.pi*(freq_start+i*k/2.0)*i + initial_phase_rad) for i in t] # just for testing #plt.plot(t,x) return x def Spectrogram_test_routine(): """ ######################################################################### # # Testing # """ # let's take a time series x(t)which was sampled at 100 Hz f_samp = 100 # Hz # Number of data points N = 2*128 # an example time series: noise + linear chirp x = np.random.rand(N) \ +np.array(Chirp_N(15.0,30.0,f_samp,N,np.pi/3.0)) # let the time point of the first data sample be the current time delta_time = dt.timedelta(seconds=1.0/f_samp) # in seconds start_time = dt.datetime.now() times = [start_time+k*delta_time for k in range(len(x))] plt.plot(times,x) plt.title('Time domain waveform: random noise + linear chirp 15=>30 Hz') plt.ylabel('value') plt.xlabel('Time (sec)') date_formatter = mdates.DateFormatter('%S.%f') ax = plt.gca() ax.xaxis.set_major_formatter(date_formatter) # calculate the spectrogram err, x_indxs, f, ampl = Spectrogram(x, fs=f_samp, \ windowfunc='hanning', \ nperwindow=16, \ nfft=64, \ slidewindowbyn = 4, \ sgram_type='ampl') #### Plot the data and set axis labels fig, ax = plt.subplots() ## Show the image with automatic color scaling #im = ax.imshow(ampl, origin='lower',interpolation='none', aspect='auto') # Show the image with manually set limits in Z data #im = ax.imshow(ampl, origin='lower',interpolation='none', aspect='auto'), \ # vmin=1.5, vmax=3.0) im = ax.imshow(ampl, origin='lower',interpolation='none', aspect='auto') # set axis labels ax.set_title('random noise + linear chirp 15=>30 Hz') ax.set_xlabel('time (sec)') ax.set_ylabel('frequency (Hz)') #### Set the tics and tick labels on the frequency axis # There are ampl.shape[0] rows, # each row corresponds to a frequency in the array freqs # in image dimensions, y axis is the number of rows, so # the value (ampl.shape[0] - 1) in the y axis corresponds to # the frequency f[-1] # If a tick needs to be shown at a specific frequency f_spec, # tick position is f_spec/f[-1] * (ampl.shape[0] - 1) # ax.set_yticks(y_tick_positions) can be used to set the ticks # ax.set_yticklabels(list(tick_labels)) can be used to set tick labels # Set up the tick positions freq_spacing = 5 # i.e., by 5 Hz y_tick_freqs = np.arange(0,f[-1],freq_spacing) y_tick_positions = (y_tick_freqs-f[0])/(f[-1]-f[0]) * (ampl.shape[0] - 1) # Let's show the ticks... ax.set_yticks(y_tick_positions) # ... and label them with the respective entries ax.set_yticklabels(y_tick_freqs) #### Set the tics and tick labels on the x (or time) axis # There are ampl.shape[1] columns, # each column corresponds to a value on x axis (or a time point) # so that x_value_for_the_i_th_ampl = x[x_indxs[i]] # or time_point_for_the_i_th_ampl = times[x_indxs[i]] # in image dimensions, x axis is the number of columns, so # the value (ampl.shape[1] - 1) in the x axis corresponds to # x[x_indxs[-1]], or in this example: times[x_indxs[-1]] # If a tick needs to be shown at a specific x_value (or time_point), # the corresponding tick position is # ( (x_value-x[x_indxs[0]]) # /(x[x_indxs[-1]]-x[x_indxs[0]]) ) * (ampl.shape[1] - 1) # or # ( (time_point-times[x_indxs[0]]).total_seconds() # /(times[x_indxs[-1]]-times[x_indxs[0]]).total_seconds() ) # * (ampl.shape[1] - 1) # ax.set_xticks(x_tick_positions) can be used to set the ticks # ax.set_xticklabels(list(tick_labels)) can be used to set tick labels time_spacing = dt.timedelta(microseconds = 100000) x_tick_times = [ times[x_indxs[0]].replace(microsecond=0) ] while ((x_tick_times[len(x_tick_times)-1]) < times[x_indxs[-1]]): x_tick_times.append((x_tick_times[len(x_tick_times)-1] + time_spacing)) time_span = times[x_indxs[-1]] - times[x_indxs[0]] x_tick_positions = [(i - times[x_indxs[0]]).total_seconds()/time_span.total_seconds() * (ampl.shape[1] - 1) \ for i in x_tick_times] #x_tick_labels = [i.strftime("%b.%d,%H:%M") for i in x_tick_times] #x_tick_labels = [i.strftime("%b.%d,%Hh") for i in x_tick_times] x_tick_labels = [ "{0:06.3f}".format(i.second+i.microsecond/1e6) for i in x_tick_times] # Let's show the ticks... ax.set_xticks(x_tick_positions) # ... and label them with the respective entries ax.set_xticklabels(x_tick_labels, rotation=90) # Set the viewing area to display the image # set y axis limits ax.set_ylim([ 0, ampl.shape[0] - 1 ]) # set x axis limits ax.set_xlim([ 0 , ampl.shape[1] - 1 ]) #### add colorbar cbar = ax.figure.colorbar(im, ax=ax) # set the label of the colorbar cbar.ax.set_ylabel('amplitude', rotation=-90, va="bottom") #### Show the plot plt.tight_layout() plt.show() ######################################################################### #Spectrogram_test_routine()