Source code for noise_floor

#!/usr/bin/env python

#########################################################################
#
#	This Program has definitions to calculate the fourier transform and noise floor
#	of a set of data from a given file.

#	Copyright (C) 2016 Aman Abhishek Tiwari, Indian Institute of Technology, Kanpur
#
#	This program is free software: you can redistribute it and/or modify
#	it under the terms of the GNU General Public License as published by
#	the Free Software Foundation, either version 3 of the License, or
#	(at your option) any later version.
#
#	This program is distributed in the hope that it will be useful,
#	but WITHOUT ANY WARRANTY; without even the implied warranty of
#	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#	GNU General Public License for more details.
#
#	You should have received a copy of the GNU General Public License
#	along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
########################################################################## 

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import sys
from scipy.signal import argrelextrema

#########################################################################

[docs]def greater (arr, comp):
	""" Compares all the elements of the an array to the number comp and generates a list which is greater than the number comp.
	
	- **parameters**, **types**, **return** and **return types**::
		:param arr: Array on which the operation is to be done.
		:param comp: The number to which the array elements is to be compared.
		:type arr: floats[]
		:type arr: float
		:return: float[]
	"""
	result = []
	for idx in range (len (arr)):
		if arr[idx] > comp:
			result [idx] = arr[idx]
	return result

#########################################################################
#########################################################################

[docs]def noise_floor_distribution (data_points) :
	y_range = 1000
	x_range = 100
	interval = len (data_points) / x_range
	result = np.zeros ((x_range,y_range))
	for i in range (x_range):
		data = np.array (\
			data_points [int (i * interval) : int ((i + 1) * interval)])
		ymin = min (data)
		ymax = max (data)
		y_interval = (ymax - ymin) / y_range
		for j in range (y_range) :
			total_p = len (data)
			result [i][j] = (data > (j * y_interval)).sum() * 100 / total_p
	
	#for row in result:
		#for elt in row:
			#print (elt)
	
	plt.figure (2)
	sp = plt.subplot (111)
	sp.set_title ("Noise Distribution below different levels.")
	sp.set_xlabel ("No. of Points")
	sp.set_ylabel ("% of points below a dB in fourier transform.")
	#x = np.array (range(len (result[400])))
	for i in range (x_range):
		sp.plot (result[i])

	plt.grid()
	return

#########################################################################
#########################################################################

[docs]def noise_floor (data_points) :
	"""Calculates the noise floor on the given fourier transform and plots it
	on the original fourier transform generated by the fourier transform module.

	- **parameters**, **types**, **return** and **return types**::
		:param data_points: Array of which noise floor is to be calculated.
		:type arr: floats[]

	"""
	plt.figure (1)
	y_range = 1000
	x_range = 1000
	interval = len (data_points) / x_range
	#result = np.zeros ((x_range,y_range))
	results = []
	for i in range (x_range):
		data = data_points [int (i * interval) : int ((i + 1) * interval)]
		#print (type (data))
		ymin = min (data)[0]
		ymax = max (data)[0]
		y_interval = (ymax - ymin) / y_range
		total_points = len (data)
		print (i)
		for j in range (y_range) :
			#result [i][j] = (data > (j * y_interval)).sum() * 100 / total_points
			tmp = ()
			total = 0
			for elt in data:
				total += elt[0] >= (j * y_interval)
				
			
			y_tmp = j * y_interval
			if total * 100 / total_points < 20:
				for elt in data:
					if elt[0] >= y_tmp:
						results.append (elt)
				break
				#print ("within if ", results)
	
	
	#for row in result:
		#for elt in row:
			#print (elt)
	
	#x = np.array (range(len (result[400])))
	
	x = []
	y = []
	for elt in results:
		x.append (elt[1])
		y.append (elt[0])

	sp2 = plt.subplot(211)
	sp2.plot (x, y, ".", label="Noise Floor")
	legend = sp2.legend(loc='upper right', shadow=True)
	sp2.grid()
	return

#########################################################################

[docs]def fourier_transform (datafile, freq, fc = 5e6):
	"""Plots the fourier transform of the given data, calculates the noise floor
	and also plots the original data in time domain along with the data being 
	passed through a high-pass filter.
	
	- **parameters**, **types**, **return** and **return types**::
		:param datafile: File from which data has to be read. This file should contain just the data points.
		:param freq: Acquisition frequency.
		:param fc: Cutoff frequency for the high-pass filter.
		:type datafile: str
		:type freq: str
		:type fc: float
	"""

	infile = open (datafile, 'r')
	freq = float (freq)

	y = []
	for line in infile:
		y.append (float (line))

	y = np.array (y)

	y_out = np.fft.fft (y)
	y_out_log = 20 * np.log10 (np.absolute (y_out) / 8191)
	#y_out_log = np.absolute (y_out)
	#y_out_log = 20 * np.absolute (y_out) / 8191
	plt.figure (1)
	x = np.array (range(len (y))) / float (freq)
	final = []
	
	###############################

	sp2 = plt.subplot(211)
	sp2.set_title ("Fourier transform of Data.")
	sp2.set_xlabel ("Frequency (s^-1)")
	sp2.set_ylabel ("dB")
	sp2.grid()
	x_out = np.fft.fftfreq (len (x), d = 1 / freq)
	
	for idx in range (int (len (x) / 2)):
		final.append (([y_out_log[idx], x_out[idx]]))
	sp2.plot (x_out[1:len (x_out)/2], y_out_log[1:len (x_out)/2], ".", label="Fourier Transform")
	legend = sp2.legend(loc='upper right', shadow=True)
	#print (min (final)[0])
	#print (type (final))
	###############################

	n = 1.0
	H = 1 / (1 + (x_out/ fc)**(2 * n))

	yflt = np.real (np.fft.ifft (H * y_out))
	sp3 = plt.subplot(212)
	sp3.set_title ("Data plot")
	sp3.set_xlabel ("Time (s)")
	sp3.set_ylabel ("Value")
	sp3.plot (x, y, label="Original Data")
	sp3.plot (x, yflt, label=("Filtered Data" + str(fc)))
	legend = sp3.legend(loc='upper right', shadow=True)
	print (max (yflt));
	print (min (yflt));
	noise_floor_distribution (y_out_log[:len (x_out)/2])
	noise_floor (final)
	sp3.grid()
	plt.show()
	return

if __name__ == '__main__':
    fourier_transform()
#########################################################################