From freelancer, 2 Weeks ago, written in Python.
Embed
  1. #!/usr/bin/env python3
  2.  
  3. import numpy as np
  4. import matplotlib.pyplot as plt
  5.  
  6.  
  7. '''signal definitions'''
  8. sample_rate = 200.0     # sample rate in samples/second
  9. signal_frequency = 5.0  # signal frequency in hertz
  10. signal_amplitude = 1.0  # signal amplitude in volt
  11. noise_stddev = 0.5      # noise standard deviation in volt
  12.  
  13.  
  14. '''generate time vector 1 second long'''
  15. sample_interval = 1 / sample_rate
  16. time = np.arange(0.0, 1.0, sample_interval)
  17. samples_number = len(time)
  18.  
  19.  
  20. '''generate one clean and one noisy (mean free AWGN noise) sine wave signal'''
  21. signal = signal_amplitude * np.sin(2 * np.pi * signal_frequency * time)
  22. noise = np.random.normal(0.0, noise_stddev, samples_number)
  23. signal_noisy = signal + noise
  24.  
  25.  
  26. '''calculate number of frequency bins for one-sided spectrum'''
  27. bin_number = int(np.ceil(samples_number / 2))
  28. frequency_bins = np.arange(bin_number)
  29.  
  30.  
  31. '''calculate FFT and normalize for one-sided spectrum'''
  32. signal_spectrum = np.fft.fft(signal_noisy) / bin_number
  33. signal_spectrum = signal_spectrum[range(bin_number)]
  34.  
  35.  
  36. '''generate graphs, plot absolute of complex spectrum'''
  37. figure, axis = plt.subplots(3, 1)
  38.  
  39. axis[0].plot(time, signal)
  40. axis[0].set_xlabel('Time (seconds)')
  41. axis[0].set_ylabel('Amplitude (volt)')
  42.  
  43. axis[1].plot(time, signal_noisy)
  44. axis[1].set_xlabel('Time (seconds)')
  45. axis[1].set_ylabel('Amplitude (volt)')
  46.  
  47. axis[2].plot(frequency_bins, abs(signal_spectrum), 'r')
  48. axis[2].set_xlabel('Frequency (hertz)')
  49. axis[2].set_ylabel('|Y(freq)|')
  50.  
  51. plt.show()