perception_scale/scale_plot.py at main · LXP-Never/perception_scale · GitHub

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
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
# -*- coding:utf-8 -*-
# Author:凌逆战 | Never
# Date: 2022/5/25
"""

"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter


def hz2mel(hz):
    """ Hz to Mels """
    return 2595 * np.log10(1 + hz / 700.0)


def hz2erb_percepNet(hz):
    """ Hz to Erbs """
    return 9.265 * np.log(1 + hz / (24.7 * 9.265))


def hz2erb_1983(f):
    """ 中心频率f(Hz) f to ERB(Hz) """
    f = f / 1000.0
    return (6.23 * (f ** 2)) + (93.39 * f) + 28.52


def hz2erb_1990(f):
    """ 中心频率f(Hz) f to ERB(Hz) """
    return 24.7 * (4.37 * f / 1000 + 1.0)


def hz2erb_1998(f):
    """ 中心频率f(Hz) f to ERB(Hz)
        hz2erb_1990 和 hz2erb_1990_2 的曲线几乎一模一样
        M. Slaney, Auditory Toolbox, Version 2, Technical Report No: 1998-010, Internal Research Corporation, 1998
        http://cobweb.ecn.purdue.edu/~malcolm/interval/1998-010/
    """
    return 24.7 + (f / 9.26449)


def Hz2erb_2002(f):
    """ [Hohmann2002] Equation 16 """
    EarQ = 9.265  # _ERB_Q
    minBW = 24.7  # minBW
    return EarQ * np.log(1 + f / (minBW * EarQ))


def Hz2erb_matlab(f):
    """ Convert Hz to ERB number """
    n_erb = 11.17268 * np.log(1 + (46.06538 * f) / (f + 14678.49))
    return n_erb


def Hz2erb_other(f):
    """ Convert Hz to ERB number """
    n_erb = 21.4 * np.log10(1 + 0.00437 * f)
    return n_erb


# freq2erb_matlab 和 Hz2erb_2002相差无几
# erb_1990和erb_1998相差无几

if __name__ == "__main__":
    fs = 16000
    hz = np.linspace(0, fs // 2, fs // 2)
    erb_1983 = hz2erb_1983(hz)
    erb_1990 = hz2erb_1990(hz)
    erb_1998 = hz2erb_1998(hz)
    erb_2002 = Hz2erb_2002(hz)
    Hz2erb_matlab = Hz2erb_matlab(hz)
    Hz2erb_other = Hz2erb_other(hz)
    mel = hz2mel(hz)
    hz2erb_percepNet = hz2erb_percepNet(hz)

    # plt.plot(hz, erb_1983, label="erb_1983")
    # plt.plot(hz, erb_1990, label="erb_1990")
    # plt.plot(hz, erb_1998, label="erb_1998")
    plt.plot(hz, erb_2002, label="erb_2002")
    plt.plot(hz, Hz2erb_matlab, label="Hz2erb_matlab")
    plt.plot(hz, Hz2erb_other, label="Hz2erb_other")
    # plt.plot(hz, mel, label="mel")
    plt.plot(hz, hz2erb_percepNet, label="hz2erb_percepNet")

    plt.legend()  # 显示图例
    plt.xlabel("Center frequency (Hz)", fontsize=12)  # x轴的名字
    plt.ylabel("Equivalent Rectangular Bandwidth (Hz)", fontsize=12)

    plt.xticks(fontsize=10)  # x轴的刻度
    plt.yticks(fontsize=10)

    plt.xlim(0, fs // 2)  # 坐标轴的范围
    plt.ylim(0)


    def formatnum(x, pos):
        return '$%.1f$' % (x / 1000)


    formatter = FuncFormatter(formatnum)
    # plt.gca().xaxis.set_major_formatter(formatter)
    # plt.gca().yaxis.set_major_formatter(formatter)
    plt.grid(linestyle='--')
    plt.tight_layout()
    plt.show()