#!/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)');