forked from jaakkopasanen/AutoEq
-
Notifications
You must be signed in to change notification settings - Fork 0
/
check_min_phase.py
61 lines (51 loc) · 1.85 KB
/
check_min_phase.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
# -*- coding: utf-8 -*-
import os
from glob import glob
import numpy as np
import matplotlib.pyplot as plt
from frequency_response import FrequencyResponse
def fft(x, fs):
nfft = len(x)
df = fs / nfft
f = np.arange(0, fs - df, df)
X = np.fft.fft(x)
X_mag = 20 * np.log10(np.abs(X))
return f[0:int(np.ceil(nfft/2))], X_mag[0:int(np.ceil(nfft/2))]
def main():
fs = 48000
f_res = 60
input_dir = os.path.join('oratory1990', 'data', 'onear')
glob_files = glob(os.path.join(input_dir, '**', '*.csv'), recursive=True)
for input_file_path in glob_files:
fr = FrequencyResponse.read_from_csv(input_file_path)
fr.equalization = fr.raw
fr.raw = np.array([])
mp = fr.minimum_phase_impulse_response(fs=fs, f_res=f_res, normalize=False)
mp = np.concatenate([mp, np.zeros(fs//10 - len(mp))])
f_mp, mp = fft(mp, fs)
f_mp[0] = 0.1
mp = FrequencyResponse(name='Minimum phase', frequency=f_mp, raw=mp)
mp.center()
lp = fr.linear_phase_impulse_response(fs=fs, f_res=f_res, normalize=False)
lp = np.concatenate([lp, np.zeros(fs//10 - len(lp))])
f_lp, lp = fft(lp, fs)
f_lp[0] = 0.1
lp = FrequencyResponse(name='Linear phase', frequency=f_lp, raw=lp)
lp.center()
fig, ax = plt.subplots()
fig.set_size_inches(15, 10)
plt.plot(fr.frequency, fr.equalization)
plt.plot(mp.frequency, mp.raw, '.-')
plt.plot(lp.frequency, lp.raw, '.-')
plt.legend(['Raw', 'Minimum phase', 'Linear phase'])
plt.semilogx()
plt.xlabel('Frequency (Hz)')
plt.xlim([1, 20000])
plt.ylabel('Gain (dBr)')
plt.ylim([-20, 20])
plt.title(fr.name)
plt.grid(True, which='major')
plt.grid(True, which='minor')
plt.show()
if __name__ == '__main__':
main()