| 1 | #=======================================================================
|
| 2 | # Direct and inverse FFT examples (including frequency analysis,
|
| 3 | # spectrogram generation and "brutal" band pass filtering)
|
| 4 | #
|
| 5 | # Cesare Brizio, 3 January 2024
|
| 6 | #
|
| 7 | # Based on several examples, no creative merit on my part except
|
| 8 | # reconciling a few variable / array names among the examples.
|
| 9 | #
|
| 10 | #=======================================================================
|
| 11 | import numpy as np
|
| 12 | from matplotlib import pyplot as plt
|
| 13 | from scipy.io import wavfile
|
| 14 | from scipy.fft import fft, fftfreq
|
| 15 | from scipy.fft import rfft, rfftfreq
|
| 16 | from scipy.fft import irfft
|
| 17 | from scipy import signal
|
| 18 |
|
| 19 |
|
| 20 | SAMPLE_RATE = 44100 #Samples per second
|
| 21 | DURATION = 5 #Seconds
|
| 22 | FREQUENCY_TEST = 2 #Hertz
|
| 23 |
|
| 24 | def generate_sine_wave(freq, sample_rate, duration):
|
| 25 | x = np.linspace(0, duration, sample_rate * duration, endpoint=False)
|
| 26 | frequencies = x * freq
|
| 27 | #2pi because np.sin takes radians
|
| 28 | y = np.sin((2 * np.pi) * frequencies)
|
| 29 | return x, y
|
| 30 |
|
| 31 |
|
| 32 | #======================================================
|
| 33 | #======================================================
|
| 34 | #
|
| 35 | # GENERATE AND PLOT A SINE WAVE
|
| 36 | #
|
| 37 | #======================================================
|
| 38 | #======================================================
|
| 39 |
|
| 40 | #Generate a 2 hertz sine wave that lasts for 5 seconds
|
| 41 | x, y = generate_sine_wave(FREQUENCY_TEST, SAMPLE_RATE, DURATION)
|
| 42 | plt.rcParams['figure.figsize'] = [12, 7]
|
| 43 | plt.plot(x, y)
|
| 44 | # displaying the title
|
| 45 | plt.title("2 Hz sine Wave generated by generate_sine_wave()\nSample rate = "+str(SAMPLE_RATE)+" Hz, Duration = "+str(DURATION)+" sec\nClose window for next step",
|
| 46 | fontsize=14,
|
| 47 | fontweight="bold",
|
| 48 | color="green")
|
| 49 | plt.ylabel('SIN value')
|
| 50 | plt.xlabel('sec')
|
| 51 | plt.show()
|
| 52 |
|
| 53 |
|
| 54 | #======================================================
|
| 55 | #======================================================
|
| 56 | #
|
| 57 | # GENERATE, NORMALIZE, PLOT AND SAVE AS A WAVEFILE
|
| 58 | # THE SUM OF TWO SINE WAVES (WILL BE BAND-STOP-FILTERED
|
| 59 | # IN SUBSEQUENT STEPS)
|
| 60 | #
|
| 61 | #======================================================
|
| 62 | #======================================================
|
| 63 | FREQUENCY_GOOD = 400
|
| 64 | FREQUENCY_NOISE = 4000
|
| 65 | _, nice_tone = generate_sine_wave(FREQUENCY_GOOD, SAMPLE_RATE, DURATION)
|
| 66 | _, noise_tone = generate_sine_wave(FREQUENCY_NOISE, SAMPLE_RATE, DURATION)
|
| 67 | noise_tone = noise_tone * 0.3
|
| 68 |
|
| 69 | mixed_tone = nice_tone + noise_tone
|
| 70 |
|
| 71 | normalized_tone = np.int16((mixed_tone / mixed_tone.max()) * 32767)
|
| 72 |
|
| 73 | plt.rcParams['figure.figsize'] = [12, 7]
|
| 74 | LIST_SLICE = 1000
|
| 75 | DURATION_LIST_SLICE = round((1/SAMPLE_RATE)*LIST_SLICE,3)
|
| 76 | FREQ_GOOD_CYCLES_PER_SLICE = round(DURATION_LIST_SLICE*FREQUENCY_GOOD,3)
|
| 77 | plt.plot(normalized_tone[:LIST_SLICE])
|
| 78 | # displaying the title
|
| 79 | plt.title("Normalized tone from the sum of 2 sine waves (will be saved as mysinewave.wav)\n\u00ABGood\u00BB = "+str(FREQUENCY_GOOD)+" Hz, \u00ABNoise\u00BB = "+str(FREQUENCY_NOISE)+" Hz, Duration = first "+str(LIST_SLICE)+" entries at a sample rate of "+str(SAMPLE_RATE)+" Hz\nClose window for next step",
|
| 80 | fontsize=14,
|
| 81 | fontweight="bold",
|
| 82 | color="green")
|
| 83 | plt.ylabel('Samples')
|
| 84 | plt.xlabel('Entries at a sample rate of '+str(SAMPLE_RATE)+' Hz - 1000 entries = '+str(DURATION_LIST_SLICE)+' sec - 1000 entries = '+str(FREQ_GOOD_CYCLES_PER_SLICE)+' cycles at '+str(FREQUENCY_GOOD)+' Hz')
|
| 85 | plt.show()
|
| 86 |
|
| 87 | #Remember SAMPLE_RATE = 44100 Hz is our playback rate
|
| 88 | wavfile.write("C:\mysinewave.wav", SAMPLE_RATE, normalized_tone)
|
| 89 |
|
| 90 |
|
| 91 | #======================================================
|
| 92 | #======================================================
|
| 93 | #
|
| 94 | # PERFORM AND DISPLAY DIRECT FFT (FROM TIME DOMAIN TO
|
| 95 | # FREQUENCY DOMAIN) OF THE NORMALIZED COMBINATION
|
| 96 | # OF TWO SINE WAVES (WILL BE BAND-STOP-FILTERED
|
| 97 | # IN SUBSEQUENT STEPS). THE FFT WILL INCLUDE THE ENTIRE
|
| 98 | # COMPLEX INPUT, INCLUDING NEGATIVE FREQUENCIES
|
| 99 | #
|
| 100 | #======================================================
|
| 101 | #======================================================
|
| 102 |
|
| 103 | #Number of samples in normalized_tone
|
| 104 | N = SAMPLE_RATE * DURATION
|
| 105 |
|
| 106 | yf = fft(normalized_tone)
|
| 107 | xf = fftfreq(N, 1 / SAMPLE_RATE)
|
| 108 |
|
| 109 | plt.rcParams['figure.figsize'] = [12, 7]
|
| 110 | plt.plot(xf, np.abs(yf))
|
| 111 | # displaying the title
|
| 112 | plt.title("Direct FFT of the normalized tone\nFull-Range ( -"+str(SAMPLE_RATE/2)+" Hz to +"+str(SAMPLE_RATE/2)+" Hz) Pressure-Frequency Analysis\nClose window for next step",
|
| 113 | fontsize=14,
|
| 114 | fontweight="bold",
|
| 115 | color="green")
|
| 116 | plt.ylabel('FFT value (sound pressure)')
|
| 117 | plt.xlabel('Frequency')
|
| 118 | plt.show()
|
| 119 |
|
| 120 |
|
| 121 |
|
| 122 | #======================================================
|
| 123 | #======================================================
|
| 124 | #
|
| 125 | # PERFORM AND DISPLAY DIRECT FFT (FROM TIME DOMAIN TO
|
| 126 | # FREQUENCY DOMAIN) OF THE NORMALIZED COMBINATION
|
| 127 | # OF TWO SINE WAVES (WILL BE BAND-STOP-FILTERED
|
| 128 | # IN SUBSEQUENT STEPS). THE FFT WILL INCLUDE ONLY
|
| 129 | # REAL INPUT (EXCLUDING NEGATIVE FREQUENCIES)
|
| 130 | #
|
| 131 | #======================================================
|
| 132 | #======================================================
|
| 133 |
|
| 134 | #Note the extra 'r' at the front
|
| 135 | yf = rfft(normalized_tone)
|
| 136 | xf = rfftfreq(N, 1 / SAMPLE_RATE)
|
| 137 | plt.rcParams['figure.figsize'] = [12, 7]
|
| 138 | plt.plot(xf, np.abs(yf))
|
| 139 | # displaying the title
|
| 140 | plt.title("Direct FFT of the normalized tone (REAL INPUT, does not compute the negative frequency terms)\nRange (0 Hz to +"+str(SAMPLE_RATE/2)+" Hz) Pressure-Frequency Analysis\nClose window for next step",
|
| 141 | fontsize=14,
|
| 142 | fontweight="bold",
|
| 143 | color="green")
|
| 144 | plt.ylabel('FFT value (sound pressure)')
|
| 145 | plt.xlabel('Frequency')
|
| 146 | plt.show()
|
| 147 |
|
| 148 |
|
| 149 | #======================================================
|
| 150 | #======================================================
|
| 151 | #
|
| 152 | # READ THE WAVEFILE AND, USING THE WELCH() FUNCTION
|
| 153 | # FROM THE SIGNAL PACKAGE, PLOT THE POWER SPECTRUM
|
| 154 | # IN dBFS.
|
| 155 | # Welch’s method computes an estimate of the power
|
| 156 | # spectral density by dividing the data into
|
| 157 | # overlapping segments, computing a modified
|
| 158 | # periodogram for each segment and averaging the
|
| 159 | # periodograms.
|
| 160 | #
|
| 161 | #======================================================
|
| 162 | #======================================================
|
| 163 |
|
| 164 | segment_size = 512
|
| 165 |
|
| 166 | fs, x = wavfile.read('C:\mysinewave.wav')
|
| 167 | x = x / 32768.0 # scale signal to [-1.0 .. 1.0]
|
| 168 |
|
| 169 | noverlap = segment_size / 2
|
| 170 | f, Pxx = signal.welch(x, # signal
|
| 171 | fs=fs, # sample rate
|
| 172 | nperseg=segment_size, # segment size
|
| 173 | window='hann', # window type to use e.g.hamming or hann
|
| 174 | nfft=segment_size, # num. of samples in FFT
|
| 175 | detrend=False, # remove DC part
|
| 176 | scaling='spectrum', # return power spectrum [V^2]
|
| 177 | noverlap=noverlap) # overlap between segments
|
| 178 |
|
| 179 | # set 0 dB to energy of sine wave with maximum amplitude
|
| 180 | ref = (1/np.sqrt(2)**2) # simply 0.5 ;)
|
| 181 | p = 10 * np.log10(Pxx/ref)
|
| 182 |
|
| 183 | plt.rcParams['figure.figsize'] = [12, 7]
|
| 184 |
|
| 185 | fill_to = -150 * (np.ones_like(p)) # anything below -150dB is irrelevant
|
| 186 | plt.fill_between(f, p, fill_to )
|
| 187 | plt.xlim([f[2], f[-1]])
|
| 188 | plt.ylim([-150, 6])
|
| 189 | # plt.xscale('log') # uncomment if you want log scale on x-axis
|
| 190 | plt.grid(True)
|
| 191 | # displaying the title
|
| 192 | plt.title("Direct FFT of the normalized tone (REAL INPUT, does not compute the negative frequency terms)\nRange (0 Hz to +"+str(SAMPLE_RATE/2)+" Hz) Pressure-Frequency Analysis\nClose window for next step",
|
| 193 | fontsize=14,
|
| 194 | fontweight="bold",
|
| 195 | color="green")
|
| 196 | plt.xlabel('Frequency, Hz')
|
| 197 | plt.ylabel('Power spectrum, dBFS')
|
| 198 | plt.show()
|
| 199 |
|
| 200 |
|
| 201 |
|
| 202 | #======================================================
|
| 203 | #======================================================
|
| 204 | #
|
| 205 | # READ THE WAVEFILE AND, USING THE SPECTROGRAM()
|
| 206 | # FUNCTION FROM THE SIGNAL PACKAGE, PLOT THE
|
| 207 | # TIME/FREQUENCY SPECTROGRAM.
|
| 208 | #
|
| 209 | #======================================================
|
| 210 | #======================================================
|
| 211 |
|
| 212 | f, t, Sxx = signal.spectrogram(x, SAMPLE_RATE)
|
| 213 | plt.pcolormesh(t, f, Sxx, shading='nearest') # shading may be flat, nearest, gouraud or auto
|
| 214 | plt.rcParams['figure.figsize'] = [12, 7]
|
| 215 | # displaying the title
|
| 216 | plt.title("Spectrogram of the normalized tone (REAL INPUT, does not compute the negative frequency terms)\nRange (0 Hz to +"+str(SAMPLE_RATE/2)+" Hz) Pressure-Frequency Analysis\nClose window for next step",
|
| 217 | fontsize=14,
|
| 218 | fontweight="bold",
|
| 219 | color="green")
|
| 220 | plt.ylabel('Frequency [Hz]')
|
| 221 | plt.xlabel('Time [sec]')
|
| 222 | plt.show()
|
| 223 |
|
| 224 |
|
| 225 | #======================================================
|
| 226 | #======================================================
|
| 227 | #
|
| 228 | # BAND-STOP FILTERING (BRUTAL, NO WINDOWING FUNCTION)
|
| 229 | #
|
| 230 | # Above, scipy.fft.rfftfreq() was used to return the
|
| 231 | # Discrete Fourier Transform sample frequencies.
|
| 232 | # The xf float array, returned by scipy.fft.rfftfreq()
|
| 233 | # contains the frequency bin centers in cycles per unit
|
| 234 | # of the sample spacing (with zero at the start). For
|
| 235 | # instance, if the sample spacing is in seconds, then
|
| 236 | # the frequency unit is cycles/second.
|
| 237 | #
|
| 238 | # The yf float array was returned by scipy.fft.rfft(),
|
| 239 | # a function that computes the 1-D n-point discrete
|
| 240 | # Fourier Transform (DFT) for real input.
|
| 241 | #
|
| 242 | # By setting to zero the yf values corresponding to
|
| 243 | # a given xf bin, we are brutally silencing the range
|
| 244 | # of frequencies subsumed by that bin.
|
| 245 | #
|
| 246 | # A more sophisticated approach to filtering, based
|
| 247 | # on windowing functions such as the one exemplified
|
| 248 | # above. A brutal surrogate of the windowing function
|
| 249 | # is setting to zero also a few bins adjacent to the
|
| 250 | # target bin.
|
| 251 | #
|
| 252 | #======================================================
|
| 253 | #======================================================
|
| 254 |
|
| 255 | #The maximum frequency is half the sample rate
|
| 256 | points_per_freq = len(xf) / (SAMPLE_RATE / 2)
|
| 257 |
|
| 258 | #Our target (noise) frequency is 4000 Hz
|
| 259 | target_idx = int(points_per_freq * 4000)
|
| 260 | #Also a small range of adjacent bins is set to zero
|
| 261 | yf[target_idx - 1 : target_idx + 2] = 0
|
| 262 |
|
| 263 | plt.rcParams['figure.figsize'] = [12, 7]
|
| 264 | plt.plot(xf, np.abs(yf))
|
| 265 | # displaying the title
|
| 266 | plt.title("Effects of the \u00ABbrutal\u00BB band stop filtering of the normalized tone \nRange (0 Hz to +"+str(SAMPLE_RATE/2)+" Hz) Pressure-Frequency Analysis\nClose window for next step",
|
| 267 | fontsize=14,
|
| 268 | fontweight="bold",
|
| 269 | color="green")
|
| 270 | plt.ylabel('FFT value (sound pressure)')
|
| 271 | plt.xlabel('Frequency')
|
| 272 | plt.show()
|
| 273 |
|
| 274 |
|
| 275 | #======================================================
|
| 276 | #======================================================
|
| 277 | #
|
| 278 | # RESTORE THE SINE WAVE (NOW NOT INCLUDING NOISE) BY
|
| 279 | # PERFORMING AN INVERSE FFT (REAL VALUES ONLY), DISPLAY
|
| 280 | # THE CLEAN WAVE AND SAVE IT AS A WAVEFILE.
|
| 281 | #
|
| 282 | #======================================================
|
| 283 | #======================================================
|
| 284 |
|
| 285 | new_sig = irfft(yf)
|
| 286 |
|
| 287 | plt.rcParams['figure.figsize'] = [12, 7]
|
| 288 | plt.plot(new_sig[:LIST_SLICE])
|
| 289 | # displaying the title
|
| 290 | plt.title("Normalized tone after the band-stop filtering of the noise wave (will be saved as clean.wav)\nDuration = first "+str(LIST_SLICE)+" entries at a sample rate of "+str(SAMPLE_RATE)+" Hz\nClose window to end program",
|
| 291 | fontsize=14,
|
| 292 | fontweight="bold",
|
| 293 | color="green")
|
| 294 | plt.ylabel('Samples')
|
| 295 | plt.xlabel('Entries at a sample rate of '+str(SAMPLE_RATE)+' Hz - 1000 entries = '+str(DURATION_LIST_SLICE)+' sec - 1000 entries = '+str(FREQ_GOOD_CYCLES_PER_SLICE)+' cycles at '+str(FREQUENCY_GOOD)+' Hz')
|
| 296 | plt.show()
|
| 297 |
|
| 298 | norm_new_sig = np.int16(new_sig * (32767 / new_sig.max()))
|
| 299 |
|
| 300 |
|
| 301 | wavfile.write("C:\clean.wav", SAMPLE_RATE, norm_new_sig)
|
