Created
January 16, 2020 19:49
-
-
Save gavv/41db26aa76737e72f70112293d7ddae4 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!env python3 | |
import numpy as np | |
import matplotlib.pyplot as plt | |
from matplotlib.pyplot import * | |
from scipy.signal import * | |
from scipy.fftpack import * | |
from optparse import OptionParser | |
import scipy.io.wavfile | |
import scipy | |
import math | |
import sys | |
def Spectrum(s): | |
Ftest = scipy.fftpack.fft( s ) | |
n = round(s.shape[0]/2) | |
xf = np.linspace(0.0, Fs/2.0, n) | |
return xf, 20*np.log10(np.abs(Ftest[0:n])/n) | |
if __name__ == "__main__": | |
parser = OptionParser() | |
parser.add_option( "-f", "--file", action="store", | |
default="/tmp/sin.wav", | |
type="string", dest="dest_file", | |
help="Where to store retrieved signal.") | |
(options, args) = parser.parse_args() | |
Fs = 48000 | |
F = 440.0 | |
SNR = 140 # dB | |
A = [8000, 1000] | |
freq=[440.0, 23e3] | |
signal_power = np.square(np.linalg.norm(A)) | |
noise_sigma = np.sqrt(signal_power / math.pow( 10, SNR/20 )) | |
def fan(time): | |
f = sum([a*np.sin(2*np.pi*f*time) for a,f in zip(A, freq)]) | |
f = f + np.random.normal(0, noise_sigma) | |
return f | |
t = np.arange( 10.0, step = 1/Fs ) | |
y = np.array([[fan(_),0] for _ in t], np.int16) | |
# y = np.zeros([t.shape[0],2],np.int16) | |
# y[round(44100*2.5), 0] = 16384 | |
scipy.io.wavfile.write(options.dest_file, Fs, y) | |
y = y.T | |
plt.figure() | |
plt.subplot(211) | |
plt.plot(*Spectrum(y[0])) | |
plt.grid() | |
plt.xlim([0, 24e3]) | |
plt.subplot(212) | |
plt.plot(y[0]) | |
plt.show() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!env python3 | |
import numpy | |
import pylab | |
import scipy.io.wavfile | |
original = scipy.io.wavfile.read('mkh10/original.wav' )[1][:,0] | |
resampler_cpy = scipy.io.wavfile.read('mkh10/resampler_cpy.wav')[1][:,0] | |
resampler_sc1 = scipy.io.wavfile.read('mkh10/resampler_sc1.wav')[1][:,0] | |
original = original / 2**15 | |
resampler_cpy = resampler_cpy / 2**31 | |
resampler_sc1 = resampler_sc1 / 2**31 | |
# resampler cuts first and last 1024 samples | |
original = original[1024:] | |
original = original[:len(original)-1024] | |
# sox cuts some more at the end | |
original = original[:len(resampler_cpy)] | |
f, [p1, p2, p3] = pylab.subplots(nrows=3, ncols=1, sharex=True, sharey=True) | |
p1.plot(original, '-', label='original') | |
p1.plot(resampler_cpy, '.', label='resampler dummy (copy)') | |
p1.plot(resampler_sc1, '-', label='resampler with scaling=1') | |
p1.legend(loc='lower left') | |
p2.plot(original - resampler_cpy, label='original - resampler dummy') | |
p2.legend(loc='lower left') | |
p3.plot(original - resampler_sc1, label='original - resampler scaling=1') | |
p3.legend(loc='lower left') | |
pylab.show() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!env python3 | |
import numpy | |
import pylab | |
import scipy.fftpack | |
import scipy.io.wavfile | |
def spectrum(s): | |
rate = 48000 | |
Ftest = scipy.fftpack.fft( s ) | |
n = round(s.shape[0]/2) | |
xf = numpy.linspace(0.0, rate/2.0, n) | |
return xf, 20*numpy.log10(numpy.abs(Ftest[0:n])/n) | |
original = scipy.io.wavfile.read('mkh10/original.wav' )[1][:,0] | |
resampler_cpy = scipy.io.wavfile.read('mkh10/resampler_cpy.wav')[1][:,0] | |
resampler_sc1 = scipy.io.wavfile.read('mkh10/resampler_sc1.wav')[1][:,0] | |
original = original / 2**15 | |
resampler_cpy = resampler_cpy / 2**31 | |
resampler_sc1 = resampler_sc1 / 2**31 | |
# resampler cuts first and last 1024 samples | |
original = original[1024:] | |
original = original[:len(original)-1024] | |
# sox cuts some more at the end | |
original = original[:len(resampler_cpy)] | |
f, [p1, p2, p3] = pylab.subplots(nrows=3, ncols=1, sharex=True, sharey=True) | |
p1.plot(*spectrum(original), label='original') | |
p1.legend(loc='upper center') | |
p2.plot(*spectrum(original), label='original') | |
p2.plot(*spectrum(resampler_cpy), label='resampler dummy (copy)') | |
p2.legend(loc='upper center') | |
p3.plot(*spectrum(original), label='original') | |
p3.plot(*spectrum(resampler_sc1), label='resampler with scaling=1') | |
p3.legend(loc='upper center') | |
pylab.show() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment