Bunch of signal processing examples by Julius Orion Smith III, converted from Matlab.

git-svn-id: https://svn.ssec.wisc.edu/repos/geoffc/Python@130 f45efd6a-04c6-4179-886d-a06c9f63a260

Applying_Blackman_Window.py

+#!/usr/bin/env python
+
+import numpy as np
+from numpy import cos,pi
+import scipy as sp
+from pylab import *
+
+
+""" Applying the Blackman taper window """
+
+M = 64
+T = 1;                       # Set sampling rate to 1
+A = 1;                       # Sinusoidal amplitude
+phi = 0;                     # Sinusoidal phase
+
+n = np.arange(M)   # Discrete time axis
+f = 0.25 + 0.5/M   # Move frequency up 1/2 bin
+
+x = cos(2.*pi*n*f*T)    # The data
+
+""" Various windows """
+
+window = np.blackman     # The window
+#window = np.hanning      # The window
+#window = np.hamming      # The window
+#window = np.bartlett     # The window
+#window = np.kaiser       # The window
+
+w = window(M)
+
+zpf = 8               # zero-padding factor
+x_Pad = concatenate((x,np.zeros((zpf-1)*M)))  # zero-padded data
+w_pad = concatenate((w,np.zeros((zpf-1)*M)))  # zerp-padded window
+
+xw = w_pad*x_Pad    # The windowed data
+
+X = np.fft.fft(xw)  # Fourier spectrum
+magX = abs(X)       # Spectral magnitude
+spec = 10.*log10(conj(X)*X) # Spectral magnitude in dB
+
+nfft = zpf*M        # Padded data length
+fni = np.arange(0.,nfft)/float(nfft-1)  # Normalized frequency axis
+fninf = fni - 0.5   # Centered frequency axis
+
+
+""" Windowed, zero-padded data """
+
+figure()
+subplot(3,1,1)
+plot(xw,'-k')
+title(window.func_name+' Windowed, Zero-Padded, Sampled Sinusoid')
+xlabel('Time (samples)')
+ylabel('Amplitude')
+text(-50,1,'a)')
+
+""" Plot spectral magnitude: """
+
+subplot(3,1,2)
+plot(fni,magX,'-k'); grid(b=True)
+title('Smoothed, Interpolated, Spectral Magnitude'); 
+xlabel('Normalized Frequency (cycles per sample))')
+ylabel('Magnitude (Linear)')
+text(-.11,40,'b)')
+
+""" Same thing on a dB scale: """
+
+subplot(3,1,3);
+plot(fninf,np.fft.fftshift(spec),'-k'); 
+xlim(-0.5,0.5)
+ylim(-40.,40.)
+grid(b=True)
+title('Smoothed, Interpolated, Spectral Magnitude (dB)'); 
+xlabel('Normalized Frequency (cycles per sample))'); 
+ylabel('Magnitude (dB)')
+text(-.6,40,'b)')
+
+show()
+
+#################################
+
+#% Windowed, zero-padded data:
+#n = [0:M-1];          % discrete time axis
+#f = 0.25 + 0.5/M;     % frequency
+#xw = [w .* cos(2*pi*n*f),zeros(1,(zpf-1)*M)]; 
+
+#% Smoothed, interpolated spectrum:
+#X = fft(xw);                                    
+
+#% Plot time data:
+#subplot(2,1,1); 
+#plot(xw);
+#title('Windowed, Zero-Padded, Sampled Sinusoid');
+#xlabel('Time (samples)'); 
+#ylabel('Amplitude');
+#text(-50,1,'a)'); 
+
+#% Plot spectral magnitude:
+#spec = 10*log10(conj(X).*X);  % Spectral magnitude in dB
+#spec = max(spec,-60*ones(1,nfft)); % clip to -60 dB
+#subplot(2,1,2);
+#plot(fninf,fftshift(spec),'-'); 
+#axis([-0.5,0.5,-60,40]); 
+#title('Smoothed, Interpolated, Spectral Magnitude (dB)'); 
+#xlabel('Normalized Frequency (cycles per sample))'); 
+#ylabel('Magnitude (dB)'); grid;
+#text(-.6,40,'b)');

FFT_Not_So_Simple_Sinusoid.py

+#!/usr/bin/env python
+
+import numpy as np
+from numpy import cos,pi
+import scipy as sp
+from pylab import *
+
+
+"""
+Example 2: = Example 1 with frequency between bins
+"""
+
+N = 64;                      # Must be a power of two
+T = 1;                       # Set sampling rate to 1
+A = 1;                       # Sinusoidal amplitude
+phi = 0;                     # Sinusoidal phase
+
+f = 0.25 + 0.5/N;   # Move frequency up 1/2 bin
+
+n = np.arange(N);            # Discrete time axis
+
+x = cos(2.*pi*n*f*T); # Signal to analyze
+X = np.fft.fft(x);          # Spectrum
+
+""" Plot time data: """
+
+figure();
+subplot(3,1,1);
+plot(n,x,'*k');        
+title('Sinusoid tuned NEAR 1/4 the Sampling Rate')
+xlabel('Time (samples)'); 
+ylabel('Amplitude');
+text(-8,1,'a)');
+
+# Interpolated sinusoid for comparison
+ni = np.arange(0.,N-1,0.1)     # Interpolated time axis
+plot(ni,A*cos(2*pi*ni*f*T+phi),'-k')
+
+""" Plot spectral magnitude: """
+
+magX = abs(X);
+fn = np.arange(0.,N)/float(N)  # Normalized frequency axis
+subplot(3,1,2);
+stem(fn,magX,linefmt='k-', markerfmt='ko'); grid(b=True) 
+xlabel('Normalized Frequency (cycles per sample))'); 
+ylabel('Magnitude (Linear)');
+text(-.11,40,'b)');
+
+""" Same thing on a dB scale: """
+
+spec = 20.*log10(magX); # Spectral magnitude in dB
+subplot(3,1,3);
+plot(fn,spec,'--ok');grid(b=True)
+xlim(0.,1.)
+ylim(-10.,30.)
+xlabel('Normalized Frequency (cycles per sample))'); 
+ylabel('Magnitude (dB)');
+text(-.11,50,'c)');
+
+""" Plot the periodic extension of the time-domain signal """
+
+#figure()
+#plot(concatenate((x,x)),'--ok');
+#title('Time Waveform Repeated Once');
+#xlabel('Time (samples)'); ylabel('Amplitude');
+
+show()

FFT_Simple_Sinusoid.py

+#!/usr/bin/env python
+
+import numpy as np
+from numpy import cos,pi
+import scipy as sp
+from pylab import *
+
+"""
+Example 1: FFT of a DFT-sinusoid
+"""
+
+""" Parameters: """
+
+N = 64;                      # Must be a power of two
+T = 1;                       # Set sampling rate to 1
+A = 1;                       # Sinusoidal amplitude
+phi = 0;                     # Sinusoidal phase
+f = 0.25;                    # Frequency (cycles/sample)
+n = np.arange(N);            # Discrete time axis
+x = A*cos(2*pi*n*f*T+phi);   # Sampled sinusoid
+X = np.fft.fft(x);           # Spectrum
+
+""" Plot time data: """
+
+figure();
+subplot(3,1,1);
+plot(n,x,'*k');        
+title('Sinusoid at 1/4 the Sampling Rate')
+xlabel('Time (samples)'); 
+ylabel('Amplitude');
+text(-8,1,'a)');
+
+# Interpolated sinusoid for comparison
+ni = np.arange(0.,N-1,0.1)     # Interpolated time axis
+plot(ni,A*cos(2*pi*ni*f*T+phi),'-k')
+
+""" Plot spectral magnitude: """
+
+magX = abs(X);
+fn = np.arange(0.,N)/float(N)  # Normalized frequency axis
+subplot(3,1,2);
+stem(fn,magX,linefmt='k-', markerfmt='ko'); grid(b=True) 
+xlabel('Normalized Frequency (cycles per sample))'); 
+ylabel('Magnitude (Linear)');
+text(-.11,40,'b)');
+
+""" Same thing on a dB scale: """
+
+spec = 20.*log10(magX); # Spectral magnitude in dB
+subplot(3,1,3);
+plot(fn,spec,'--ok');grid(b=True)
+xlim(0.,1.)
+ylim(-350.,50.)
+xlabel('Normalized Frequency (cycles per sample))'); 
+ylabel('Magnitude (dB)');
+text(-.11,50,'c)');
+
+""" Plot the periodic extension of the time-domain signal """
+
+#figure()
+#plot(concatenate((x,x)),'--ok');
+#title('Time Waveform Repeated Once');
+#xlabel('Time (samples)'); ylabel('Amplitude');
+
+show()

FFT_Zero_Padded_Sinusoid.py

+#!/usr/bin/env python
+
+import numpy as np
+from numpy import cos,pi
+import scipy as sp
+from pylab import *
+
+
+"""
+Example 3: = Example 2 with zero padding
+"""
+
+N = 64;                      # Must be a power of two
+T = 1;                       # Set sampling rate to 1
+A = 1;                       # Sinusoidal amplitude
+phi = 0;                     # Sinusoidal phase
+
+f = 0.25 + 0.5/N   # Move frequency up 1/2 bin
+
+n = np.arange(N)   # Discrete time axis
+
+zpf = 8            # zero-padding factor
+
+"""
+Original signal will form 1/zpf of padded signal,
+with zero padding making up the rest.
+"""
+x = concatenate((cos(2.*pi*n*f*T),np.zeros((zpf-1)*N)))  # zero-padded
+n = np.arange(zpf*N)   # Discrete time axis for padded data
+
+X = np.fft.fft(x)  # Spectrum
+magX = abs(X)      # magnitude spectrum
+
+""" Plot time data: """
+
+figure()
+subplot(3,1,1)
+plot(n,x,'-k')
+title('Zero-padded Sampled Sinusoid')
+xlabel('Time (samples)')
+ylabel('Amplitude')
+text(-8,1,'a)')
+
+""" Plot spectral magnitude: """
+
+subplot(3,1,2)
+nfft = zpf*N       # FFT size = new frequency grid size
+fni = np.arange(0.,nfft)/float(nfft-1)  # Normalized frequency axis
+# with interpolation, we can use solid lines '-':
+plot(fni,magX,'-k'); grid(b=True)
+xlabel('Normalized Frequency (cycles per sample))')
+ylabel('Magnitude (Linear)')
+text(-.11,40,'b)')
+
+""" Same thing on a dB scale: """
+
+subplot(3,1,3)
+spec = 20.*log10(magX) # Spectral magnitude in dB
+#spec = max(spec,-40*np.ones(size(spec))) # clip below at -40 dB:
+plot(fni,spec,'-k');grid(b=True)
+xlim(0.,1.)
+ylim(-40.,40.)
+xlabel('Normalized Frequency (cycles per sample))')
+ylabel('Magnitude (dB)')
+text(-.11,50,'c)')
+
+show()

Hann_Window_Spectrum_Analysis.py

+#!/usr/bin/env python
+
+import numpy as np
+from numpy import cos,pi
+import scipy as sp
+from pylab import *
+
+""" Hann-Window: Spectal Analysis """
+
+# Analysis parameters:
+M = 33         # Window length 
+N = 64         # FFT length (zero padding factor near 2)
+
+# Signal parameters:
+wxT = 2.*pi/4.  # Sinusoid frequency (rad/sample)
+A = 1.          # Sinusoid amplitude
+phix = 0.       # Sinusoid phase
+
+# Compute the signal x:
+n = np.arange(N)    # time indices for sinusoid and FFT
+x = A * exp(1j*wxT*n+phix) # complex sine [1,j,-1,-j...]
+
+# Compute Hann window:
+nm = np.arange(M)   # time indices for window computation
+
+# Hann window = "raised cosine", normalization (1/M)
+# chosen to give spectral peak magnitude at 1/2:
+w = (1./M) * (cos((pi/M)*(nm-(M-1.)/2.)))**2.  
+
+wzp = concatenate((w,np.zeros(N-M))) # zero-pad out to the length of x
+xw = x * wzp          # apply the window w to signal x
+
+# Compute the spectrum and its alternative forms:
+Xw = np.fft.fft(xw);              # FFT of windowed data
+fn = np.arange(N)/float(N)  # Normalized frequency axis
+spec = 20.*log10(abs(Xw))  # Spectral magnitude in dB
+phs = angle(Xw)           # Spectral phase in radians
+phsu = unwrap(phs)        # Unwrapped spectral phase
+
+# Compute heavily interpolated versions for comparison:
+Nzp = 16                   # Zero-padding factor
+Nfft = N*Nzp               # Increased FFT size
+xwi = concatenate((xw,np.zeros(Nfft-N))) # New zero-padded FFT buffer
+Xwi = np.fft.fft(xwi)             # Compute interpolated spectrum
+fni = np.arange(Nfft)/float(Nfft)  # Normalized frequency axis
+speci = 20.*log10(abs(Xwi)); # Interpolated spec mag (dB)
+phsi = angle(Xwi);          # Interpolated phase
+phsiu = unwrap(phsi);       # Unwrapped interpolated phase
+
+figure();
+subplot(2,1,1);
+
+plot(fn,abs(Xw),'*k') 
+plot(fni,abs(Xwi),'-k')
+xlim(0.,1.)
+#ylim(-100.,0.)
+title('Spectral Magnitude')
+xlabel('Normalized Frequency (cycles per sample))')
+ylabel('Amplitude (linear)')
+
+subplot(2,1,2);
+
+# Same thing on a dB scale
+plot(fn,spec,'*k')
+plot(fni,speci,'-k')
+xlim(0.,1.)
+ylim(-100.,0.)
+title('Spectral Magnitude (dB)'); 
+xlabel('Normalized Frequency (cycles per sample))'); 
+ylabel('Magnitude (dB)');
+
+#cmd = ['print -deps ', 'specmag.eps']; disp(cmd); eval(cmd);
+#disp 'pausing for RETURN (check the plot). . .'; pause
+
+figure();
+subplot(2,1,1);
+plot(fn,phs,'*k')
+plot(fni,phsi,'-k')
+xlim(0.,1.)
+ylim(-4.,4.)
+title('Spectral Phase')
+xlabel('Normalized Frequency (cycles per sample))')
+ylabel('Phase (rad)')
+grid(b=True)
+
+subplot(2,1,2)
+plot(fn,phsu,'*k')
+plot(fni,phsiu,'-k')
+xlim(0.,1.)
+ylim(-14.,0.)
+title('Unwrapped Spectral Phase')
+xlabel('Normalized Frequency (cycles per sample))')
+ylabel('Phase (rad)')
+grid(b=True)
+
+show()

Hann_Windowed_Complex_Sinusoid.py

+#!/usr/bin/env python
+
+import numpy as np
+from numpy import cos,pi
+import scipy as sp
+from pylab import *
+
+
+""" Hann-Windowed Complex Sinusoid """
+ 
+# Analysis parameters:
+M = 31         # Window length 
+N = 64         # FFT length (zero padding factor near 2)
+
+# Signal parameters:
+wxT = 2.*pi/4.  # Sinusoid frequency (rad/sample)
+A = 1.          # Sinusoid amplitude
+phix = 0.       # Sinusoid phase
+
+# Compute the signal x:
+n = np.arange(N)    # time indices for sinusoid and FFT
+x = A * exp(1j*wxT*n+phix) # complex sine [1,j,-1,-j...]
+
+# Compute Hann window:
+nm = np.arange(M)   # time indices for window computation
+
+# Hann window = "raised cosine", normalization (1/M)
+# chosen to give spectral peak magnitude at 1/2:
+w = (1./M) * (cos((pi/M)*(nm-(M-1.)/2.)))**2.  
+
+wzp = concatenate((w,np.zeros(N-M))) # zero-pad out to the length of x
+xw = x * wzp          # apply the window w to signal x
+
+figure()
+subplot(1,1,1)
+
+# Display real part of windowed signal and Hann window
+plot(n,wzp,'-k')
+plot(n,real(xw),'*b')
+title(['Hann Window and Windowed, Zero-Padded, Sinusoid (Real Part)'])
+xlabel('Time (samples)')
+ylabel('Amplitude')
+
+show()

Use_Blackman_Window.py

+#!/usr/bin/env python
+
+import numpy as np
+from numpy import cos,pi
+import scipy as sp
+from pylab import *
+
+
+""" Define the Blackman taper window """
+
+M = 64
+w = np.blackman(M)
+
+figure()
+subplot(3,1,1)
+plot(w,'*')
+xlim(0.,M)
+ylim(0.,1.)
+title('Blackman Window')
+xlabel('Time (samples)')
+ylabel('Amplitude')
+text(-8,1,'a)')
+
+""" Also show the window transform """
+
+subplot(3,1,2)
+zpf = 8                       # zero-padding factor
+xw = concatenate((w,np.zeros((zpf-1)*M)))  # zero-padded window
+Xw = np.fft.fft(xw)           # Blackman window transform
+spec = 20.*log10(abs(Xw))     # Spectral magnitude in dB
+spec = spec - max(spec)       # Normalize to 0 db max
+nfft = zpf*M                  # Padded data length
+
+fni = np.arange(0.,nfft)/float(nfft)  # Normalized frequency axis
+#plot(fni,spec,'-k')
+plot(spec,'-k')
+#xlim(0.,1.)
+ylim(-100.,10.)
+grid(b=True)
+xlabel('Normalized Frequency (cycles per sample))'); 
+ylabel('Magnitude (dB)')
+text(-.12,20,'b)'); 
+
+""" Replot interpreting upper bin numbers as frequencies<0 """
+
+subplot(3,1,3);
+#specnf = np.fft.fftshift(spec)
+specnf = np.roll(spec,nfft/2)
+fninf = fni - 0.5
+plot(fninf,specnf,'-k')
+#plot(specnf,'-k')
+xlim(-0.5,0.5)
+ylim(-100.,10.)
+grid(b=True)
+xlabel('Normalized Frequency (cycles per sample))'); 
+ylabel('Magnitude (dB)');
+text(-.62,20,'c)'); 
+
+
+
+show()

correlTests.py

+#!/usr/bin/env python
+
+import numpy as np
+from numpy import cos,pi
+#import scipy as sp
+from scipy import signal
+from pylab import *
+import time
+
+
+N = 256
+M = 256
+fullSize = N+M-1
+sameSize = max(M, N)
+validSize = max(M, N) - min(M, N) + 1
+
+nReals = 10000
+x = np.random.randn(nReals,N)
+x_vec = np.zeros(N)
+
+# Using scipy FFT convolution
+
+"""
+     mode -- a flag indicating the size of the output
+              'valid'  (0): The output consists only of those elements that
+                              are computed by scaling the larger array with all
+                              the values of the smaller array.
+              'same'   (1): The output is the same size as the largest input
+                              centered with respect to the 'full' output.
+              'full'   (2): The output is the full discrete linear convolution
+                              of the inputs. (Default)
+"""
+
+#c_1 = np.zeros(fullSize)
+#c_vec = np.zeros(fullSize)
+
+#cumTime=0.
+#for reals in range(nReals):
+    #x_vec = x[reals]
+    #t1 = time.clock()
+    #c_vec = signal.fftconvolve(x_vec,x_vec[::-1],mode='full')
+    #t2=time.clock()
+    #cumTime += t2-t1
+    #c_1 += c_vec
+
+#print "Time to perform FFT correlation: ",cumTime
+#c_1 /= N
+
+# Using regular convolution
+
+"""
+    mode : {'full', 'valid', 'same'}, optional
+        'full':
+          By default, mode is 'full'.  This returns the convolution
+          at each point of overlap, with an output shape of (N+M-1,). At
+          the end-points of the convolution, the signals do not overlap
+          completely, and boundary effects may be seen.
+
+        'same':
+          Mode `same` returns output of length ``max(M, N)``.  Boundary
+          effects are still visible.
+
+        'valid':
+          Mode `valid` returns output of length
+          ``max(M, N) - min(M, N) + 1``.  The convolution product is only given
+          for points where the signals overlap completely.  Values outside
+          the signal boundary have no effect.
+"""
+
+#c_2 = np.zeros(fullSize)
+#c_vec = np.zeros(fullSize)
+
+#cumTime=0.
+#for reals in range(nReals):
+    #x_vec = x[reals]
+    #t1 = time.clock()
+    #c_vec = np.convolve(x_vec,x_vec[::-1],mode='full')
+    #t2=time.clock()
+    #cumTime += t2-t1
+    #c_2 += c_vec
+
+#print "Time to perform regular correlation: ",cumTime
+#c_2 /= N
+
+# Using manual FFT convolution
+
+c_3 = np.zeros(fullSize)
+c_vec = np.zeros(fullSize)
+
+cumTime=0.
+x_shape = np.array(x[0].shape)
+size = x_shape+x_shape-1
+
+for reals in range(nReals):
+    x_vec = x[reals]
+    t1 = time.clock()
+    X_vec = np.fft.fftn(x_vec,size)
+    X_vec *= conj(X_vec)
+    t2=time.clock()
+    cumTime += t2-t1
+    c_3 += X_vec
+
+figure()
+plot(np.fft.fftshift(X_vec))
+
+c_3 /= N
+c_3 = np.fft.ifftn(c_3,size)
+c_3 /= c_3[0]
+c_3 = np.fft.fftshift(c_3)
+print "Time to perform manual FFT correlation: ",cumTime
+
+#######
+
+#plot(c_1)
+#plot(c_2)
+figure()
+plot(c_3)
+
+show()
+
+del(x)