Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import numpy as np
- import matplotlib.pyplot as plt
- from scipy import signal
- from scipy.io import wavfile
- import librosa
- import librosa.display
- from scipy.constants import pi
- import scipy
- gamma = 0.3
- k_max = 15 # k_max+1 kerrosta
- # read in the audio file
- fs,x = wavfile.read('/Users/mac/OneDrive - TUNI.fi/SGN/Audio/Projekti/rhythm_birdland.wav')
- x = x[0:fs]
- f, t, F = signal.stft(x, fs, nperseg=1000)
- amp = 2 * np.sqrt(2)
- # print(np.zeros((2,3)))=[0 0 0;0 0 0]
- fig1 = plt.figure(1)
- plt.pcolormesh(t, f, np.abs(F), vmin=0, vmax=amp*10)
- plt.title('STFT Magnitude')
- plt.ylabel('Frequency [Hz]')
- plt.xlabel('Time [sec]')
- plt.colorbar()
- plt.show()
- W = (np.abs(F))**(2*gamma)
- H = 0.5*W
- P = H
- alpha = 0.3 # (np.var(P))/(np.var(H)+np.var(P)) # 0.3
- delta = np.zeros((len(f-2),len(t-2))) # Delta tulee yhden jaljessa -> deltaa on k_max kerrosta
- print(" ")
- print(np.size(x), fs)
- print(P.shape)
- print(delta.shape)
- print(" ")
- for k in range(0,k_max):
- H_im1 = np.c_[ np.zeros(len(f)), H ]
- H_ip1 = np.c_[ H, np.zeros(len(f))]
- print(H_im1)
- P_hm1 = np.vstack([ np.transpose(np.zeros(len(t))), P ])
- P_hp1 = np.vstack([ P, np.transpose(np.zeros(len(t))) ])
- print(P_hm1)
- print(H_im1.shape)
- print(H_ip1.shape)
- print(P_hm1.shape)
- print(P_hp1.shape)
- delta = alpha*(H_im1-2*H+H_ip1)/4 - (1-alpha)*(P_hm1-2*P+P_hp1)/4
- if k == k_max:
- H_km1 = H
- P_km1 = P
- P = W - H
- H = min(max(H + delta,0), W)
- if H_km1 < P_km1:
- H = 0
- P = W
- else:
- P = 0
- H = W
- print(" ")
- print(H.shape)
- print(P.shape)
- print(" ")
- h = scipy.signal.istft(H**(1/(2*gamma))*np.exp(1.j*np.angle(F)))
- print("Eka ISTFT tehty. Koko h = ", np.size(h))
- p = scipy.signal.istft(P**(1/(2*gamma))*np.exp(1.j*np.angle(F)))
- print("Tonen ISTFT tehty. Koko p = ", np.size(p))
- # s = audio["samples"]
- # s_q = np.square(s)
- # y = h+p
- # SNR_h = 10*np.log*()
- fig2 = plt.figure(2)
- plt.pcolormesh(t, f, np.abs(H-P), vmin=0, vmax=amp)
- plt.title('STFT Magnitude')
- plt.ylabel('Frequency [Hz]')
- plt.xlabel('Time [sec]')
- plt.colorbar()
- fig3 = plt.figure(3)
- plt.pcolormesh(t, f, np.abs(P), vmin=0, vmax=amp)
- plt.title('STFT Magnitude')
- plt.ylabel('Frequency [Hz]')
- plt.xlabel('Time [sec]')
- plt.colorbar()
- plt.show()
- h = h[0::3]
- print("Eka ISTFT tehty. Koko h = ", h)
- print(max(h))
- # Record the audio file
- # audio_name = 'projekti_h'
- # wavfile.write(audio_name, fs, h / np.max(np.abs(h)))
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement