from __future__ import division from scipy import signal def db20(W,Nfft=None): 'Given DFT, return power level in dB or take DFT if need be' assert np.isscalar(W) or W.ndim==1 if Nfft is None: # assume W is DFT return 20*log10(abs(W)) else: # assume time-domain passed, so need DFT DFT= fft.fft(array(W).flatten(),Nfft)/sqrt(Nfft) return 20*log10(abs(DFT.flatten())) def undb20(x): 'Invert the db20. Recover amplitude' return 10**(x/20.) fig, ax = subplots() fig.set_size_inches((6,3)) Ns= 16 Nf = 256*2 freqs = arange(Nf)*2*pi/Nf w = signal.hanning(Ns,False) W = db20(w,Nf) ax.plot(freqs,W,'-b',ms=4.) ax.set_ylim(ymin = -60) ax.set_xlim(xmax = pi*1.01) ax.set_xlabel('Discrete Frequency',fontsize=14) ax.set_ylabel(r'$20\log_{10}|W|$',fontsize=18) ax.grid() ax.set_title('Hanning Window',fontsize=18) ax.annotate('',fontsize=28, xy=(76/Nf*2*pi,W[0]), xytext=(76/Nf*2*pi,W[0]-32), arrowprops={'facecolor':'b','arrowstyle':'<->'}, ) ax.text( 0.4,0.5,'Peak sidelobe level', fontsize=18, transform=ax.transAxes, bbox={'fc':'y','alpha':.3}) # fig.savefig('figure_00@.png', bbox_inches='tight', dpi=300) def peak_sidelobe(w,N=256,return_index=False, return_all=False): '''Given window function, return peak sidelobe level and bin index of all (return_all=True) or some sidelobe peaks if desired (return_index=True). Note that this method fails when the window function has no roots on the unit circle (e.g. exponential window). The return index is in units of DFT-bin (k/N). ''' assert (len(w)<=N) # need longer DFT otherwise r=np.roots(w) # find complex roots of window function r = r[np.where(np.round(abs(r),3)==1)] # keep only those on unit circle (approx) y=log(r).imag/2./pi*N # get k^th bin index y=y[y>0].astype(np.int32) # keep positive half only as integer roundoff y=np.unique(y) # dump repeated y.sort() # sort in-place W = 20*log10(abs(fft.fft(w,N))) #compute DFT # loop through slices and pick out max() as peak for that slice's sidelobe sidelobe_levels = [] sidelobe_idx =[] for s in [slice(i,j) for i,j in zip(y[:-1],y[1:])]: imx= s.start+W[s].argmax() # bin index of max peak= W[imx]-W[0] # relative to global peak sidelobe_levels.append( peak ) # store sidelobe level for later sidelobe_idx.append(imx/N) # ... with corresponding bin if return_all: return zip(sidelobe_levels, sidelobe_idx) if return_index: return (sidelobe_levels[0], sidelobe_idx[0]) return sidelobe_levels[0] def dftmatrix(N=32,Ns=None): 'construct DFT matrix of size N give Ns time-samples' k= np.arange(N) if Ns is None: Ns = N n = arange(Ns) U = matrix(exp(1j* 2*pi/N *k*n[:,None])) # use numpy broadcasting to create matrix return U/sqrt(N) Ns = 64 N= 512 U=dftmatrix(N=N,Ns=Ns) offset=8 # place DFT near middle of plot for readability u=array(U[:,offset]).flatten()*sqrt(N) # phase shifts w = signal.hanning(Ns,False) level,idx = peak_sidelobe(w,N,return_index=True) x0 = u*ones(Ns) x1=u*exp(1j*2*pi*arange(Ns)*(idx)) # signal on peak of sidelobe fig,axs = subplots(2,1,sharex=True,sharey=True) fig.set_size_inches((9,4)) ax=axs[0] #ax.plot(abs(fft.fft(w*(x0),N))) ax.plot(db20(w*x0,N)) ax.arrow( offset+idx*N,-60,0,60, length_includes_head=True,lw=2., head_length=14,head_width=4,fc='g',alpha=0.3) #ax.arrow( idx*N,0,0,3,length_includes_head=True,lw=1.5,head_width=2,fc='g') ax.arrow(offset,-60,0,60, length_includes_head=True, lw=2.,head_length=14,head_width=4,fc='b') #ax.legend(loc=0) ax.set_xlim(xmax=N/4.,xmin=-3) ax.set_ylim(ymax = 17,ymin=-60) ax.set_ylabel(r'$20\log_{10}|W_k|$',fontsize=18) ax.text(0.4,.5,'''The signal (green) on the first sidelobe is reduced by 31 dB during the convolution but still contaminates the other signal (blue) by the residual amount''',va='center',fontsize=12,transform=ax.transAxes); ax=axs[1] ax.plot(db20(w*x1,N)) ax.arrow( offset+idx*N,-60,0,60,length_includes_head=True,lw=2.,head_length=14,head_width=4,fc='g') ax.arrow( offset,-60,0,60,length_includes_head=True,lw=2.,head_length=14,head_width=4,fc='b',alpha=0.3) #ax.legend(loc=0) ax.set_xlim(xmax=N/4.,xmin=-3) ax.set_ylim(ymax = 17,ymin=-60) ax.set_ylabel(r'$20\log_{10}|W_k|$',fontsize=18) ax.set_xlabel('k',fontsize=16) ax.text(0.4,.6,'''In this case, the window's mainlobe peak reaches the green signal and it's the blue signal's turn to sit on the sidelobe and contaminate the green signal.''',va='center',fontsize=12,transform=ax.transAxes); # fig.savefig('figure_00@.png', bbox_inches='tight', dpi=300) fig,ax = subplots() fig.set_size_inches((7,3)) N=512 w=signal.windows.hamming(Ns) W=db20(w,N) m =10 p=np.polyfit(arange(m)/N*Ns,W[:m]-W[0]+3.01,2) # fit quadratic polynomial width = np.roots(p)[0]*2 # 3-dB beamwidth ax.plot(arange(N)/N*Ns,W-W[0]) # normalize to peak ax.set_ylim(ymin=-10) ax.set_xlim(xmax = 2) ax.vlines(width/2,0,-60,lw=2.,linestyle='--',color='g') ax.set_ylabel(r'$20\log_{10}|X|$',fontsize=22) ax.set_title(r'$ BW_{3dB}$=%3.2f bins'%width) ax.set_xlabel(r'$\frac{N_s}{N} k$',fontsize=22) ax.annotate('',fontsize=28,xy=(0,-3), xytext=(width/2,-3), arrowprops=dict(arrowstyle="<->",lw=3)) ax.annotate('',fontsize=28,xy=(1.2,0), xytext=(1.2,-3), arrowprops=dict(arrowstyle="<->",lw=3)) ax.hlines(-3,width/2,2,linestyle='--',color='g',lw=2.) ax.text( width/2/4,-5,r'$\frac{BW_{3dB}}{2}$',fontsize=22) ax.text( 1.3,-2,'-3 dB',fontsize=18) # fig.savefig('figure_00@.png', bbox_inches='tight', dpi=300) fig,ax = subplots() fig.set_size_inches((7,3)) N=256 Ns = 32 w=signal.windows.triang(Ns) W=db20(w,N) W0 = db20(exp(1j*2*pi/Ns*arange(Ns)*1/2.)*w,N)-W[0] W=W-W[0] # rescale for plot ax.plot(W,label='window') ax.plot(W0,label='half-bin shifted') scalloping_loss = W[0]-W0[0] ax.axis(ymin=-30,ymax=1,xmax=24) ax.set_title('Triangular window scalloping Loss=%3.2f'%(scalloping_loss)) ax.set_xlabel('k',fontsize=16) ax.set_ylabel(r'$20\log_{10}|W|$',fontsize=22) ax.hlines(-scalloping_loss,0,24,color='red',linestyle='--',lw=2.) ax.legend(loc=0) ax.grid() # fig.savefig('figure_00@.png', bbox_inches='tight', dpi=300)