2021-9-30 背景噪声的研究

（一）背景噪声的研究

声音如何产生？
声音是由物体振动产生声波，通过介质（空气或固体、液体）传播，并能被人的听觉器官所感知的波动现象。
最初发出振动的物体叫声源。声音以波的形式振动传播，所以声音是一种波。
信号处理的意义？
声音的本质是压力脉动，所以无论是仿真还是实验，直接测量得到的都只是每个时刻的压力值。尽管我们也可以通过时域上的压力 信号的大小定性的判断压力脉动的强弱，但是压力的时间信号太过杂乱无章，我们通常需要把时域的信号转换为频域，以便更清晰的辨识信号的特征。
时域与频域？
时域：描述信号随时间变化的坐标系。
频域：描述信号在频率特性时用到的一种坐标系，其横坐标就是频率。
对于特定几何尺度的物体，气流在一秒钟扫过它的次数就是频率，所以速度越快，意味着产生的信号频率越高；而频率与几何尺度则往往成反比的关系，即尺寸越小，产生的频率越高。相对于时域，频域的信号更容易和几何结构以及流场建立相互对应的关系。
如何转换？
傅里叶原理：任何连续测量的时域信号，都可以表示为不同频率的正弦波信号的无限叠加。傅里叶变换，则是把这些正弦波信号复原并换算到频域空间的积分变换方法。

对于离散时间序号信号$x(j),j=1,2,3,\dots ,N$,其傅里叶变换及其逆变换可以表示为：
$$Y(K)=\sum _{j=1}^{N}x(j)W_N^{(j-1)(k-1)}$$
$$x(j)=\frac{1}{n}\sum _{k=1}^{N}Y(k)W_N^{-(j-1)(k-1)}$$
其中，$W_N=e^{(-2\pi i)/N}$。
离散傅里叶变换也是从时域到频域的转换，只不过它是把一组复数$x(j)$，转换成了另一组复数$Y(K)$，就可以愉快的实现对离散数据的变换。
离散傅里叶变换好是好，只不过如果采样的数据量大了，计算量也会急剧增加。1965年，Cooley和Tukey提出了离散傅里叶变换（DFT）的快速算法，即快速傅里叶变换。
频谱
我们知道频谱的横坐标是频率，纵坐标表示该频率下的振幅，对于声音信号来说也就是压力脉动的强弱。
在matlab中实现傅里叶变换

Ts = 1/50;
t = 0:Ts:10-Ts;
x = sin(2*pi*15*t) + sin(2*pi*20*t);
plot(t,x)% 横坐标为时间，纵坐标为值
xlabel('Time(seconds)')
ylabel('Amplitude')
y = fft(x);
fs = 1/Ts;
f = (0:length(y)-1)*fs/length(y);%设置傅里叶变换y为纵轴，f为横轴
plot(f,abs(y))%最后生成频域图
xlabel('Frequency(Hz)')
ylabel('Magnitude')
title('Magnitude')
%fftshift 函数对变换执行以零为中心的循环平移
n = length(x);
fshift = (-n/2:n/2-1)*(fs/n);
yshift = fftshift(y);
plot(fshift,abs(yshift))
xlabel('Frequency (Hz)')
ylabel('Magnitude')
rng('default')
xnoise = x + 2.5*randn(size(t));%在x中混入高斯噪声
%绘画含噪声的图像
ynoise = fft(xnoise);
ynoiseshift = fftshift(ynoise);
power = abs(ynoiseshift).^2/n;
plot(fshift,power)
title('Power')
xlabel('Frequency (Hz)')
ylabel('Power')
%计算快速傅里叶
whaleFile = 'bluewhale.au';
[x,fs] = audioread(whaleFile);
whaleMoan = x(2.45e4:3.10e4);
t = 10*(0:1/fs:(length(whaleMoan)-1)/fs);

plot(t,whaleMoan)
xlabel('Time (seconds)')
ylabel('Amplitude')
xlim([0 t(end)])
m = length(whaleMoan); 
n = pow2(nextpow2(m));
y = fft(whaleMoan,n); 
f = (0:n-1)*(fs/n)/10; % frequency vector
power = abs(y).^2/n;   % power spectrum      
plot(f(1:floor(n/2)),power(1:floor(n/2)))
xlabel('Frequency (Hz)')
ylabel('Power')

达芬方程(Duffing equation)

数学表达式：
$\ddot{x}+\delta \dot{x}+\alpha x+\beta x^3=\gamma cos(wt)$
$x_1=\dot{x}$
$x_2=\ddot{x}$
我们可以变化得到：
$x_1=\dot{x}$
$x_2=-\delta \dot{x}-\alpha x-\beta x^3+\gamma cos(wt)$

matlab实现
1.写一个函数

function y = duffineq(t,x)

global delta alpha beta gamma omega

x1 = x(2);
x2 = -delta * x(2) - alpha * x(1) - beta * x(1)^3 + gamma * cos(omega * t);

y = [x1;x2];

end

2.画图

global delta alpha beta gamma omega

delta = 1;%阻尼系数,(Ns/m or kg/s)
alpha = 1;%线性刚度,(N/m)
beta = 1;%非线性刚度,(N/m^3)
gamma = 1;%驱动力的大小,(N)
omega = 1;%驱动力的角速度,(rad/s)

fsize = 18;%字号大小

x0 = [0.1 0]; %初始状态 [初始位移，初始速度]
tspan = 0:0.001:100; %仿真时长

[t, y] = ode45(@duffineq, tspan, x0);

figure
plot(y(:,1),y(:,2)) %第一列y(:,1)是x的位移,第二列y(:,2)是x的速度。
xlabel('$\it x$','Interpreter','latex','FontSize',fsize)
ylabel('$\dot{x}$', 'Interpreter','latex','FontSize',fsize)

3.运行结果

Python实现
1.代码

import numpy as np
import math
from scipy.integrate import odeint
import matplotlib.pyplot as plt


# function
def model(x, t):
    delta = 1
    alpha = 1
    beta = 1
    gamma = 1
    omega = 1
    x1 = x[1]
    x2 = -delta * x[1] - alpha * x[0] - beta * x[0] ** 3 + gamma * math.cos(omega * t)
    return [x1, x2]


# initial condition
x0 = [0.1, 0]

# time
t = np.linspace(0, 100, 100000)  # 起始点，终止点，个数。

# solve ode
y = odeint(model, x0, t)

# plot
plt.plot(y[:, 0], y[:, 1])
plt.xlabel('x')
plt.ylabel(r"$\. x $")  # $\. x $ 为x头上一点
plt.title(r"x vs. $\. x $")
plt.show()

2.运行结果