text: """
  Genetic Algorithms for Digital Signal Processing
  Created on Mon Oct 05 20:01:05 2020
  Last Edited on  Mon Oct 12 2020 by Luke Trenberth
  """
  import numpy as np
  from scipy import signal
  from scipy.fftpack import fft
  import matplotlib.pyplot as plt
  import DSP_GA as ga
  
  class DSP_Signal():
      def __init__(self, filename, fs=1024, N_Coeff=400):
          file = open(filename, "r")
          self.y_0 = []
          for line in file:
              words = line.split("  ")
              for word in words:
                  if word != "":
                      self.y_0.append(float(word))
          self.fs = fs
          self.N = len(self.y_0)
          self.N_2 = int(self.N/2)
          self.t = [x/self.fs for x in list(range(0, self.N))]
          self.f = [x*self.fs/self.N for x in list(range(0, self.N_2))]
          self.P_0 = np.var(self.y_0)
          self.FFT_0 = fft(self.y_0)
          
          self.N_Coeff = N_Coeff  # Number of coefficients
          
      #Window Filtering method for the data class
      def WF(self, GA_data):
          #GA Data: [noise_f_1, noise_f_2, width]
          
          # Window Filtering
          self.width_WF = 8 # Width of stop band, Hz
          self.band_1 = [GA_data[0] -GA_data[2]/2, GA_data[0]+GA_data[2]/2] # Define band 1 bounds
          self.band_2 = [GA_data[1] -GA_data[2]/2, GA_data[1]+GA_data[2]/2] # Define band 2 bounds
          
          self.filter1_WF = signal.firwin(self.N_Coeff+1, self.band_1, window='hann', pass_zero='bandstop', fs=self.fs) # Filter for noise frequency 1
          self.filter2_WF = signal.firwin(self.N_Coeff+1, self.band_2, window='hann', pass_zero='bandstop', fs=self.fs) # Filter for noise frequency 2
          self.filter_WF = signal.convolve(self.filter1_WF, self.filter2_WF)  # Combined filter for noise frequencies
          self.y_WF = signal.lfilter(self.filter_WF, 1, self.y_0) # Apply noise filters to original data
          self.f_WF, self.h_WF = signal.freqz(self.filter_WF, 1, fs=self.fs) #
          self.FFT_WF = fft(self.y_WF)
          return self.SNR(self.y_WF)
          
      #Parks McLellan Filtering Method
      def PM(self, GA_data, TW =3, BW=5):
          # Filter Bands for filtering frequency 1 & 2
          f_1 = GA_data[0]
          f_2 = GA_data[1]
          if len(GA_data) > 2:
              TW = GA_data[2]
          if len(GA_data) > 3:
              BW = GA_data[3]
          
          band1_PM = [0, f_1 -BW/2-TW, f_1 -BW/2, f_1+BW/2, f_1+BW/2+TW, self.fs/2]
          band2_PM = [0, f_2 -BW/2-TW, f_2 -BW/2, f_2+BW/2, f_2+BW/2+TW, self.fs/2]
          gain_PM = [1, 0, 1]
          
          # Create filters for filtering frequency 1 & 2
          filter1_PM = signal.remez(self.N_Coeff+1, band1_PM, gain_PM, fs=self.fs) # Filter frequency 1
          filter2_PM = signal.remez(self.N_Coeff+1, band2_PM, gain_PM, fs=self.fs) # Filter frequency 2
          filter_PM = signal.convolve(filter1_PM, filter2_PM) # Combined Filter
          
          self.y_PM = signal.lfilter(filter_PM, 1, self.y_0) # Filter original data in time domain
          self.f_PM, self.h_PM = signal.freqz(filter_PM, 1, fs=self.fs) # Return filter frequency response
          self.FFT_PM = fft(self.y_PM) # Filtered data frequency domain response
          return self.SNR(self.y_PM)
          
  
      #Returns a Signal to Noise Ratio for a given input Power
      def SNR (self, y):
          return self.P_0  - np.var(y)
      
      
      # Plots a Fast Fourier Transform for simple graphing
      def FFTplot(self, f, FFT, title="ECG Signal Frequency Spectrum"):
          plt.figure()
          plt.plot(f, abs(FFT)[:self.N_2])
          plt.xlabel("Frequency (Hz)")
          plt.ylabel("Voltage (uV)")
          plt.title(title)
          plt.show()
          
          
  
  
  
  
  #The GA_filter function filters an input waveform 
  def GA_filter(waveform, input_num, solutions_per_population, mating_parent_number, num_generations):
      
  
      # Defining the population size.
      pop_size = (solutions_per_population,input_num) # The population will have sol_per_pop chromosome where each chromosome has num_weights genes.
      #Creating the initial population.
      new_population = ga.create_population(pop_size)
      
      
      best_outputs = []
      for generation in range(num_generations):
          # Measuring the fitness of each chromosome in the population.
          fitness = ga.cal_pop_fitness(waveform, new_population)
          # The best result in the current iteration.
          best_outputs.append(np.max(fitness))
          # Selecting the best parents in the population for mating.
          parents = ga.select_mating_pool(new_population, fitness, 
                                            mating_parent_number)
          # Generating next generation using crossover.
          offspring_crossover = ga.crossover(parents, offspring_size=(pop_size[0]-parents.shape[0], input_num))
          # Adding some variations to the offspring using mutation.
          offspring_mutation = ga.mutation(offspring_crossover, num_mutations=2)
          # Creating the new population based on the parents and offspring.
          new_population[0:parents.shape[0], :] = parents
          new_population[parents.shape[0]:, :] = offspring_mutation
          
  #        if (generation < 20):
  #            print("{}\n {}\n\n".format(new_population, pop_fitness))
          if (generation%10 == 0 and generation != 0):
              print("{} Generations Completed".format(generation))
              
          
      # Getting the best solution after iterating finishing all generations.
      #At first, the fitness is calculated for each solution in the final generation.
      fitness = ga.cal_pop_fitness(waveform, new_population)
      # Then return the index of that solution corresponding to the best fitness.
      best_match_idx = np.where(fitness == np.max(fitness))[0]
      return new_population[best_match_idx, :], fitness[best_match_idx][0][0], best_outputs
      
  
  # Implementation of a Parks-McLellan Filter using Genetic Algorithms
  def main():
      waveform = DSP_Signal("Signal_files/enel420_grp_15.txt")
      
      # Fixed Parameters, found by trial and error s
      f_count = 2
      mating_parent_number = 3
      pop_size = 10
      num_generations = 1000
      
      # Conduct a Genetic Algorithm approximation
      best_soln, best_soln_fitness, best_outputs = GA_filter(waveform, 
                                                             f_count, pop_size, 
                                                             mating_parent_number, num_generations)    
      print("Best solution : \n", best_soln)
      print("Best solution fitness : \n", best_soln_fitness)
      plt.figure()
      plt.plot(best_outputs)
      plt.title("Fitness of ECG Signal using GA Algoritm")
      
      plt.xlabel("Number of Iterations")
      plt.ylabel("Fitness (Signal to Noise Ratio)")
      plt.grid()
      plt.show()
  
      waveform.FFTplot(waveform.f, waveform.FFT_0, title="Before filtering")
      waveform.PM(best_soln[0])
      waveform.FFTplot(waveform.f, waveform.FFT_PM, title="After Filtering")
      
  main()
  
  