# -*- coding: utf-8 -*-
"""
Inicialização para o 1º trabalho de laboratório de Sinais e Sistemas

Created on Wed Feb 26 18:38:20 2014

@author: lba
"""

def u(t):
	"""Unit step function.
	"""
	return float_(t>=0)          # float_ is the numpy standard function for conversion to float (equivalent to float64)

def delta(t):
	y = zeros(shape(t))
	y[argmin(abs(t))] = samplingrate
	return y

def tplot(x, mode='extended'):
	"""Plot a continuous-time signal with the correct time scale (symmetrical around the origin).
	If mode == 'extended' (the default), zeros are appended at both ends.
	If mode != 'extended', only the signal x is plotted, with no appended zeros.
	"""
	sizex = size(x)
	if mode == 'extended':
		x = append(zeros((sizex-1)/10), x)		# Append zeros in the beginning...
		x = append(x, zeros((sizex-1)/10))		# ...and in the end
	duration = size(x) / samplingrate
	t = timevar(duration)						# Time variable local to this function
	plot(t,x)

def tplotu(x,):
	"""Plot a continuous-time signal with the correct time scale (non-negative).
	No zeros are appended to the signal.
	"""
	sizex = size(x)
	duration = size(x) / samplingrate
	t = timevaru(duration)						# Time variable local to this function
	plot(t,x)

def fplot(xf):
	"""Plot the argument, xf, which is assumed to be in the frequency domain.
	"""
	if not all(imag(xf) == 0):
		print('fplot: The input data are not real. Cannot plot.')
		return
	sizexf = size(xf)
	f = arange(-sizexf/2, sizexf/2-.9999, 1) * 2 * pi * samplingrate / sizexf
	plot(f,xf)
	
def dplot(x, mode='extended'):
	"""Plot the argument as a discrete-time signal.
	If mode == 'extended' (the default), zeros are appended at both ends.
	If mode != 'extended', only the signal x is plotted, with no appended zeros.
	"""
	nsamples = len(x)
	if mode == 'extended':
		x = append(zeros(nsamples/10), x)
		x = append(x, zeros(nsamples/10))
		nsamples = len(x)
	n = arange(nsamples) - nsamples / 2	# Discrete time variable
	stem(n, x)

def timevar(duration):
	"""Produce a time variable in [-duration/2,duration/2) with the sampling rate defined in 'samplingrate'.
	The value duration/2 is not produced, so as to return an even, highly composite number of samples
	if samplingrate is highly composite and duration is an even integer.
	"""
	tlim = duration / 2
	t = arange(-round(tlim * samplingrate), round(tlim * samplingrate), 1) / samplingrate    # This gives exact values for integer values of t
	return t
	
def timevaru(duration):
	"""Produce a time variable in [0,duration) with the sampling rate defined in 'samplingrate'.
	The value duration is not produced, so as to return an even, highly composite number of samples
	if samplingrate is highly composite and duration is an integer.
	"""
	t = arange(0, round(duration * samplingrate), 1) / samplingrate    # This gives exact values for integer values of t
	return t
	
def play(x):
	"""Play an audio signal. The amplitude is clipped to [-1,1].
	"""
	x = clip(x, -1, 1)									# If |x|>1 this would give a strong distortion
	wavwrite('play.wav', samplingrate, (32767*x).astype(int16))
	PlaySound('play.wav',0)								# Set the second argument to 1 to returm immediately to the calling program
#	os.system('play.wav')
	os.remove('play.wav')

def seqsin(*arg):
	"""Generate a sequence of sinusoids with specified frequencies and durations.
	The arguments come in pairs: frequency and duration of each sinusoid.
	A frequency of zero corresponds to a period of silence.
	"""
	if mod(len(arg),2) != 0:
		print('seqsin: Number of arguments must be even')
		return

	t1 = .001 * samplingrate							# Duration, in samples, of the initial taper
	t2 = .01 * samplingrate							# Duration, in samples, of the final taper
	y = zeros(0)									# Initialize output array
	pos = 0										# Initialize position in which we are in the output array
	for nseg in range(int(len(arg) / 2)):
		freq = arg[2*nseg]
		duration = arg[2*nseg+1] * samplingrate			# Duration of this segment, in samples
		i = arange(t1)
		y = append(y, .5 * (1 + cos(pi / t1 * (i - t1))) * sin(2 * pi * freq / samplingrate * i))
		i = arange(t1, duration-t2)
		y = append(y, sin(2 * pi * freq / samplingrate * i))
		i = arange(duration-t2, duration)
		y = append(y, .5 * (1 + cos(pi / t2 * (i - duration + t2 + 1))) * sin(2 * pi * freq / samplingrate * i))
		pos += duration
	return y

def FourierTransform(x):
	"""Compute the Fourier transform of the signal in x.
	We use the FFT, and then adjust the output amplitude to correspond to the Fourier transform,
	taking into account the value of samplingrate.
	
	If x is not real, an error message is printed.
	
	If the size of x is not highly factorizable, this function may take quite a bit of time!
	"""
	xf = fft.fftshift(fft.fft(fft.ifftshift(x))) / samplingrate
	return xf
	
 
def convolution(x1, x2, mode=''):
	"""Compute the linear convolution of two signals.
	If mode != 'full' (the default), the two signals must have the same duration, and the result is cropped to that duration.
	If mode == 'full', the two signals may be of any durations, and the result is not cropped.
	"""

#	pdb.set_trace()

	if mode != 'full' and size(x1) != size(x2):
		print("convolution: The two signals must be of the same duration, or then you must provide the third argument with a value of 'full'.\n")
		return
	
	sizex1 = size(x1)
	sizex2 = size(x2)
	sizey = sizex1 + sizex2 - 1
	
	fftsize = 2 ** ceil(log2(sizey))			# First power of 2 that is not smaller than the size of the output
	
	x1 = append(x1, zeros(fftsize - sizex1))	# Extend the input arrays to the fft size
	x2 = append(x2, zeros(fftsize - sizex2))
	
	x1f = fft.fft(x1)						# Compute the Fourier transforms of the inputs
	x2f = fft.fft(x2)
	
	y = real(fft.ifft(x1f * x2f))			# Compute the output
	
	if mode == 'full':
		y = y[0:sizey] / samplingrate		# Set the output to the correct size for the full convolution, and to the correct amplitude
	else:
		y = y[(sizex1-1)/2 : (sizex1-1)/2 + sizex1] / samplingrate		# Truncate the output to the original size and set it to the correct amplitude
	
	return y


def sample(xc, T):
	"""Sample the signal in xc with sampling period T, which must be a multiple of 1/samplingrate.
	"""
	sampleratio = max(round(T * samplingrate), 1)
	if abs(sampleratio - T * samplingrate) > 1e-10:
		T = sampleratio / samplingrate
		print('sample: T was not a multiple of 1/samplingrate; using T=%7.5f .\n' % T)
	duration = len(xc) / samplingrate
	nsamples = 2 * floor(duration / (2 * T))
	xd = zeros(nsamples)
	for n in arange(nsamples):
		xd[n] = xc[round(duration / 2 * samplingrate - nsamples / 2 * sampleratio + n * sampleratio)]
	return xd
	
def reconstruct(yd, T):
	"""Reconstruct a continuous time signal from the samples in yd, with sampling period T, which must be a multiple of 1/samplingrate.
	We implement the reconstruction filter in the frequency domain, and therefore use circular convolution, which causes wraparound
	effects in the time domain. To reduce these effects, we compute the response with the double of the final duration, and then
	truncate it to the final duration.
	"""
#	pdb.set_trace()
	sampleratio = max(round(T * samplingrate), 1)
	if abs(sampleratio - T * samplingrate) > 1e-10:
		T = sampleratio / samplingrate
		print('reconstruct: T was not a multiple of 1/samplingrate; using T=%7.5f .\n' % T)
	nsamplesd = len(yd)				# Number of samples of yd
	nsamplesc = nsamplesd * sampleratio	# Number of samples of yc
	zerod = floor(nsamplesd / 2)			# Sample that corresponds to t=0 in yd
	zeroc = floor(nsamplesc / 2)			# Sample that corresponds to t=0 in yc
	yc = zeros(nsamplesc)
	for n in arange(nsamplesd):
		yc[(n - zerod) * sampleratio + zeroc] = yd[n] * samplingrate		# Put a delta impulse for each discrete sample
	
	# Construct reconstruction filter
#	taperfraction = .001							# Fraction of the bandwidth that is tapered
	bandwidth = round(nsamplesc / (2 * sampleratio))		# Bandwidth of the reconstruction filter in samples
#	taperwidth = round(bandwidth * taperfraction)
	hf = zeros(nsamplesc)
	hf[zeroc-bandwidth : zeroc+bandwidth+1] = 1			# Passband
#	hf[zeroc+bandwidth+1-taperwidth : zeroc+bandwidth+1] = .5 * (1 + cos(linspace(0, pi, taperwidth)))	# Taper the high end of the passband
#	hf[zeroc-bandwidth : zeroc-bandwidth+taperwidth] = .5 * (1 + cos(linspace(-pi, 0, taperwidth)))		# Taper the low end of the passband
	hf = hf * T									# Adjust the gain
	
	# Filter the sequence of deltas
	yc = real(fft.ifft(fft.fft(yc) * fft.ifftshift(hf)))
	
	return yc
		
	
						

		
	
import pdb				

from numpy import *
from matplotlib.pyplot import plot, stem
from scipy.signal import square as squarewave
import os
from scipy.io.wavfile import write as wavwrite
from scipy.io.wavfile import read as wavread
from winsound import PlaySound


samplingrate = 1000

print("\n\nSinais e Sistemas - 3º trabalho de laboratório - Controlo dum motor: inicialização concluída.\n")