当前位置：首页 > news >正文

python 处理阻尼正弦

news 来源：原创 2024/4/25 15:24:01

python 处理阻尼正弦

阻尼正弦
阻尼正弦波（英語：damped sine wave）是振幅會隨時間增長而趨向零的正弦波函數。當諧振子消耗的能量比供應的能量多，其波形即為阻尼正弦波，此函數常用在科學及工程中。
阻尼正弦可以定义为:x(t)=sin(2πft)e^(−αt)
以下代码生成正弦信号和衰减形状。通过将两个信号逐个样本相乘来生成阻尼正弦信号。

fs=8000
t=np.arange(0,1,1/fs) # crete time vector of 1 second length
f = 2.5               # 2.5 Hz 
alpha = 1             # daming factor (play around with it and see what happens)

sin  = np.sin(2*np.pi*f*t)
damp = np.exp(-alpha*t)
x = sin * damp

plt.subplot(2,1,1)
plt.plot(t,sin)
plt.ylabel('sin($t$)')

plt.subplot(2,1,2)
plt.plot(t,x, label='sin$(2 \pi \cdot ' + str(f) + ' \mathrm{Hz} \cdot t) \cdot \mathrm{e}^{- a t}$')
plt.plot(t,damp, '--', label='$\mathrm{e}^{- a t}$')
plt.legend()
plt.xlabel('$t$ in seconds')

exp：高等数学里以自然常数e为底的指数函数
numpy.exp()：返回e的幂次方，e是一个常数为2.71828

why sin = np.sin(2np.pif*t) inside the bracket it still need multiply np?

首先没有阻尼。

concert_pitch_1sec = get_sine_wave(440, 1, 8000)
half_concert_pitch_1sec = get_sine_wave(220, 1, 8000)

alternation_1sec = np.concatenate((half_concert_pitch_1sec, concert_pitch_1sec,
                                     half_concert_pitch_1sec, concert_pitch_1sec))

#plot the tune
plt.plot(alternation_1sec) # you can't see much here, it will get clear together with the next plot

ipd.Audio(alternation_1sec, rate=8000)

why inside the np.concatenate() function, half_concert_pitch_1sec, concert_pitch_1sec are repeated for twice.
why plot them display like a rectangle

下面我们可以观察到阻尼对音乐音调的影响——阻尼合成的曲调听起来更“自然”。 （此技术用于 MIDI 乐器合成）。

alternation_damped = np.concatenate((half_concert_pitch_1sec* damp, concert_pitch_1sec* damp,
                                     half_concert_pitch_1sec* damp, concert_pitch_1sec* damp))

#plot the tune
plt.plot(alternation_damped)

ipd.Audio(alternation_damped, rate=8000)

通过使用 get_sine_wave() 函数生成不同长度和频率的正弦波，您现在可以播放简单的歌曲。正如实验表中已经描述的那样，每个音符都有一个基本的 [频率]

#create musical notes 
g = get_sine_wave(196.00)
a = get_sine_wave(220)
b = get_sine_wave(246.94)
d = get_sine_wave(293.66)
c = get_sine_wave(261.63)
e = get_sine_wave(329.63)
#append together to form tune + chord 
tune = [a,d,a,g,a,b,d,(g+b+d)]
tune = np.concatenate(tune)
#plot the tune
plt.plot(tune)
#audio player for the tune 
ipd.Audio(tune,rate=8000) 

#我们可以再次观察到阻尼对音乐音调的影响。
#create dampened musical notes 
g = get_sine_wave(196.00)*damp
a = get_sine_wave(220)*damp
b = get_sine_wave(246.94)*damp
d = get_sine_wave(293.66)*damp
c = get_sine_wave(261.63)*damp
e = get_sine_wave(329.63)*damp

#create tune 
tune = [b,d,a,g,a,b,d,(g+b+d)]
tune = np.concatenate(tune)
#plot dampened tune waveform
plt.plot(tune)
#audio player for dampened tune
ipd.Audio(tune,rate=8000)

why g a b c d e have different value in get_sine_wave() but they print same plots

音符的长度由减半给出：一个二分音符的长度是一个全音符的一半，一个四分音符是半个二分音符（或一个全音符的四分之一）等。音符时长通常也可以更复杂，例如为 3/4 或 3/2。
您应该为整个音符选择合适的 fs_Hz 和持续时间。特别要确保将 fs_Hz 与十六分音符的持续时间相乘得到一个整数结果，不需要四舍五入。
我们需要一个简单的附加功能来生成更高级的歌曲是生成静音的功能。以下函数 get_silence() 生成所需长度（以秒为单位）的 0 值数组。

right_hand_notes = [
    ("E5", 1 / 16),
    ("D#5", 1 / 16),
    ("E5", 1 / 16),
    ("D#5", 1 / 16),
    ("E5", 1 / 16),
    ("B4", 1 / 16),
    ("D5", 1 / 16),
    ("C5", 1 / 16),
    ("A4", 1 / 8),
    ("Pause", 1 / 16),
    ("C4", 1 / 16),
    ("E4", 1 / 16),
    ("A4", 1 / 16),
    ("B4", 1 / 8),
    ("Pause", 1 / 16),
    ("E4", 1 / 16),
    ("G#4", 1 / 16),
    ("B4", 1 / 16),
    ("C4", 1 / 8),
    ("Pause", 1 / 16),
    ("E4", 1 / 16),
    ("E5", 1 / 16),
    ("D#5", 1 / 16),
    ("E5", 1 / 16),
    ("D#5", 1 / 16),
    ("E5", 1 / 16),
    ("B4", 1 / 16),
    ("D5", 1 / 16),
    ("C5", 1 / 16),
    ("A4", 1 / 8),
    ("Pause", 1 / 16),
    ("C4", 1 / 16),
    ("E4", 1 / 16),
    ("A4", 1 / 16),
    ("B4", 1 / 8),
    ("Pause", 1 / 16),
    ("E4", 1 / 16),
    ("C5", 1 / 16),
    ("B4", 1 / 16),
    ("A4", 1 / 4),
]

left_hand_notes = [
    ("Pause", 1 / 8),
    ("Pause", 3 / 8),
    ("A2", 1 / 16),
    ("E3", 1 / 16),
    ("A3", 1 / 16),
    ("Pause", 3 / 16),
    ("E2", 1 / 16),
    ("E3", 1 / 16),
    ("G#3", 1 / 16),
    ("Pause", 3 / 16),
    ("A2", 1 / 16),
    ("E3", 1 / 16),
    ("B3", 1 / 16),
    ("Pause", 3 / 16),
    ("Pause", 3 / 8),
    ("A2", 1 / 16),
    ("E3", 1 / 16),
    ("A3", 1 / 16),
    ("Pause", 3 / 16),
    ("E2", 1 / 16),
    ("E3", 1 / 16),
    ("G#3", 1 / 16),
    ("Pause", 3 / 16),
    ("A2", 1 / 16),
    ("E3", 1 / 16),
    ("B3", 1 / 16),
    ("Pause", 1 / 16),
]

def get_silence(length_s, sample_rate_hz):
    """Return silence for the given length at the given sample rate."""
    return np.zeros(int(length_s * sample_rate_hz))

np.zeros() 可以创建一维数组,例如np.zeros(5) -> [0. 0. 0. 0. 0.]
np.zeros() 也可以创建多维数组,例如np.zeros((5,2))
多维数组结果是:
[[0. 0.]
[0. 0.]
[0. 0.]
[0. 0.]
[0. 0.]]

现在我们可以定义一个函数来创建所需的音符或暂停，由上面对 left_hand_notes 和 right_hand_notes 的描述定义。

请注意，所需信号频率由下式计算
f=440⋅2^((note−position/12)+(octave−4))

def create_tone(note, duration):
    
    tempo = 5

    if note == "Pause":
        return get_silence(length_s=duration*tempo, sample_rate_hz=8000)

    note_position = {
        "C": -9,
        "C#": -8,
        "D": -7,
        "D#": -6,
        "E": -5,
        "F": -4,
        "F#": -3,
        "G": -2,
        "G#": -1,
        "A": 0,
        "A#": 1,
        "B": 2,
    }

    octave = int(note[-1])
    key = note[:-1]

    frequency_hz = 440 * 2 ** ((note_position[key] / 12) + (octave - 4))

    return get_sine_wave(frequency_Hz=frequency_hz, length_s=duration*tempo, fs_Hz=8000)

创建一个测试音以查看上述功能是否按预期工作。

test_tone=create_tone("A4", 1/8)

plt.plot(test_tone[1:500])
plt.title('test tone A4 (first 500 samples)')
plt.xlabel('sample index $k$')
ipd.Audio(test_tone, rate=8000)

why change the value of test_tone, the graph is still same with original

现在让我们来听一首完整的歌曲。我们可以将歌曲的双手相加得到完整的歌曲。如果您愿意，现在可以更改代码以包含阻尼。

right_hand = np.concatenate(
    [create_tone(note, duration) for note, duration in right_hand_notes]
)

left_hand = np.concatenate(
    [create_tone(note, duration) for note, duration in left_hand_notes]
)

song = left_hand + right_hand
ipd.Audio(song, rate=8000)

以下（已知）命令可用于将 WAVE 文件下载到硬盘和内存中。
请注意，如果要将变量 file_name 传递给 curl 命令，则必须将其放在大括号中。这允许您重新使用下载命令

file_name = 'music_44k.wav'
!curl https://staffwww.dcs.shef.ac.uk/people/S.Goetze/sound/{file_name} -o {file_name}

import soundfile as sf # for loading wave files from hard disk

# download
file_name = 'music_44k.wav'
print('Downloading file: ' + file_name)
!curl https://staffwww.dcs.shef.ac.uk/people/S.Goetze/sound/{file_name} -o {file_name}

# load WAVE files
mus_44k, fs44 = sf.read(file_name)

print('The sampling frequency is ' + str(fs44) + ' Hz')

# This is how our speech signal looks like
plt.plot(mus_44k);
plt.xlabel('sample index $k$');

我们使用 librosa 库来改变采样频率 $f_s$ 。以下代码首先将声音信号下采样到 $f_s=8$ kHz。通过这样做，所有高于 $f_{\mathrm{max}}=4$ kHz 的信息，即采样频率的一半，都会丢失。即使我们再次上采样到 $f_s=44.1$ kHz 的原始采样频率，丢失的信息也不会恢复。

import librosa         # we will use the library librosa here for resampling

#donwsample from 44.1 kHz to 8 kHz
mus_8k = librosa.resample(mus_44k, orig_sr=fs44, target_sr=8000); # resample to 8 kHz
#upsample from 8 kHz to 44.1 kHz
mus_44k_2 = librosa.resample(mus_8k, orig_sr=8000, target_sr=fs44); # resample to 44.1 kHz

# visualise
fig=plt.figure(figsize=(12,6)) # create a figure of size 12 x 6 inches

plt.subplot(1,3,1)
plt.specgram(mus_44k, Fs=fs44);
plt.title('original')
plt.xlabel('time $t$')
plt.ylabel('frequency $f$')

plt.colorbar(label='dB');
plt.clim(-150,0)

plt.subplot(1,3,2)
plt.specgram(mus_8k, Fs=8000);
plt.title('downsampled to $8000$ Hz')
plt.xlabel('time $t$')
plt.ylabel('frequency $f$')

plt.colorbar(label='dB');
plt.clim(-150,0)

plt.subplot(1,3,3)
plt.specgram(mus_44k_2, Fs=fs44);
plt.title('after down-/upsampling')
plt.xlabel('time $t$')
plt.ylabel('frequency $f$')

plt.colorbar(label='dB');
plt.clim(-150,0)

plt.tight_layout() # this command ensures that labelsdon't overlap

上面的三个时频表示[频谱图]可视化了下采样和上采样的效果。左侧面板显示了在 $f_s=44100 \,\mathrm{Hz}$ 处采样的原始波形文件的扇区图。您可以看到信号中存在的能量频率高达 $\frac{1}{2} f_s = 22050 \, \mathrm{Hz}$ （由 Nyquist-Shannon 采样定理： $f_{\mathrm{max}} = \frac{1}{2} f_s$ )。中间面板显示了下采样到 $f_s=8000 \,\mathrm{Hz}$ 的信号，因此频率高达 $f_{\mathrm{max}}=4000 \,\mathrm{Hz}$ 的信号能量存在于信号（参见不同的 $y$ 轴）。再次将此信号上采样到 $f_s=44100 \,\mathrm{Hz}$ 的原始采样频率（参见右图）后，您会看到频率高于 $\,\mathrm{Hz} 的信号能量$ 不会再次恢复，但会在下采样/上采样过程中丢失。
这可以通过播放音频文件来感知。