matlab/Matlab数字信号与声音处理

数字信号

函数信息与数字信号

数字信号区别与传统的函数信息的最大特点,就是其离散性,即数字信号的自变量与返回值都是离散定义的。例如,一个声音数字信号,一定有一个每秒若干次的采样,每一个时间段,声音的强度也将有一个固定精度的取值。

信号的数字化一般需要三个步骤:抽样、量化、编码

  • 一维的数字信号有时可能是波形信号,这样的信号通过拟合、插值,往往代表这一种具有光滑性与周期性的函数。
  • 一维的数字信号也有可能是声音信号,这样的信号往往会具有多种频段的信息,对应的连续函数的傅里叶变换也将在部分区段或点位出现明显更大的系数。
  • 一维的数字信号还有可能是分片光滑函数或分片常函数,也有可能是符合某种特定分布的随机函数。

信噪比

衡量误差大小的方法有很多,一种方式就是图示+观察,较 为直观可靠度也比较高。统计一致误差(最大误差)也是 一种衡量噪声强度的方法。

  • 相较来说,整体平方误差 $\Sigma[y-f(x)]^{2}$ (SE) 或平均平方误 差(MSE-mean square error ) 被更多采用。
  • 另一种基于平方误差的噪声强度的定义为信噪比 (SNRsignal to noise ratio) $: 10 \cdot \log _{10}\left(\frac{\sum f^{2}(x)}{\sum[y-f(x)]^{2}}\right)_{\circ}$ 此值与真实信号的性质有关。在相同真实信号的情况下,信噪比越高,信号质量越高。MATLAB函数调用方法snr (s, n), s​ 为信号, n为误差。

在已知噪声的分布类型与相应参数的情况下,直接以将这些噪声特性作为噪声强弱的度量,很多时候会更加可靠。 而信噪比可以作为一种衡量去噪或恢复质量的有效尺度。

例:正弦函数加噪声计算信噪比:

1
2
3
4
5
x = 0:0.001:2*pi;
f = sin(x);
f_ = f + normrnd(0, 0.1, 1, length(x));
plot(x, f, x, f_);
snr(f, f-f_) % 结果为16.6676

数字信号高斯噪声的均值滤波

因为高斯噪声是0均值的,并且噪声的幅度明显小于信号本身的强度(SNR>10即表明噪声强度不足信号强度的根号十分之一)。故采用简单的均值方法对带噪信号做光滑化处理(也称为均值滤波),就可以得到一个近似的恢复信号。

f1=smooth(y) 是一种MATLAB自带的五点均值滤波方法,返回值f1为列向量,恢复效果如下图:(SNR=22.3808)

代码为:

1
2
3
4
f1 = smooth(f_); f2 = smooth(f1); f3 = smooth(f2);
plot(x, f1, x, f2, x, f3);
legend('smooth函数的恢复结果', '2次smooth的恢复结果', '3次smooth的恢复结果')
hold off;
  • f1=conv(y,ones(1,5)*1/5,'same') 进行五点均值滤波,除边界外与smooth函数效果相同:(SNR=22.3660)
  • f2=conv(y,[1/16 1/4 3/8 1/4 1/16],'same')进行五点加权均值滤波,但本例效果一般:(SNR=20.9296)
  • f3=conv(y,ones(1,9)*1/9,'same') 进行九点均值滤波,更大的滤波范围改进了恢复效果,但滤波范围过大会导致图像过度平滑、失去应有斜率而失真:(SNR=24.3451)

恢复效果为:

傅里叶变换

连续傅里叶变换

傅里叶变换:是傅里叶级数的推广,可以从某种角度上理解成傅里叶级数的“周期”趋于无穷大的极限形式

  • 傅里叶变换定义为:$F(\omega)=\int_{-\infty}^{+\infty} f(t) \mathrm{e}^{-j \omega t} \mathrm{~d} t,$

    其中$j$为虚数单位,另外注意到, $\mathbf{e}^{-j \omega t}=\cos \omega \boldsymbol{t}-\boldsymbol{j} \cdot \sin \omega \boldsymbol{t}$

  • 傅里叶逆变换定义: $\mathrm{f}(\mathrm{t})=\frac{1}{2 \pi} \int_{-\infty}^{+\infty} F(\omega) \mathrm{e}^{j \omega t} \mathrm{~d} \omega$

  • 帕萨瓦尔等式 $\int_{-\infty}^{+\infty}|f(t)|^{2} \mathrm{~d} t=\frac{1}{2 \pi} \int_{-\infty}^{+\infty}|F(\omega)|^{2} \mathrm{~d} \omega$

与傅里叶级数类似,傅里叶变换的意义是进行时频分析, $\boldsymbol{t}$ 为时间, $\omega$ 为频率,傅里叶变换的取值往往决定了信号在什么频率上更强

Matlab函数定义

  • 函数fourier(ft, t, w) 时间域变量t,频率域变量W
  • 函数ifourier (ft, w, t)频率域变量w, 时间域变量t

离散傅里叶变换

一维离散傅里叶变换定义为$Y(k) = \sum_{i=1}^n X(i)\text{e}^{\frac{(-2\pi j)(i-1)(k-1)}{n}}, 1\leq k \leq n,$ j为虚数单位。

特别的当 $k=1$ 时, $Y(\mathbf{1})=\sum_{i=1}^{n} X(i)$ 称为离散傅里早变换的低频系数,其余系数称为高频傅里叶系数。低频系数类似于连续傅里早变换的F(0),表达了函数的总和。高频系数则在不同程度刻画函数的变化情况。

注:函数fft(x) 可以利用上式计算出向量X的离散傅里叶变换

  • 离散傅里叶变换是一种从有限长度的离散序列X(i) 到另一种有限长度离散序列Y(k) 的变换。变换的原理事实上是将序列X(i) 假想为了以n为周期的“无限长周期序列”,进而可以利用有限个特定的频率采样完成X(i)的频率表示。
  • 离散傅里叶逆变换定义为$X(i) = \frac1n\sum_{k=1}^n Y(k)\text{e}^{\frac{(2\pi j)(i-1)(k-1)}{n}}, 1\leq i \leq n,$ 这一变换与傅里叶变换完全互逆。MATLAB函数ifft(Y) 可以实现对应的功能

下图为前例正弦函数及含噪声的正弦函数的离散傅里叶变换:

代码为:

1
2
3
plot(x, abs(fft(f_)), x, abs(fft(f)));
ylim([0 20])
legend('含噪声的离散傅里叶变换', '真实信号的离散傅里叶变换')

利用离散傅里叶变换完成高斯噪声去噪

算法设计的原理即对含噪信号做离散傅里叶变换,然后仅保留大于5的系数,去除其他所有系数,再用离散傅里叶逆变换得到恢复后的数字信号,代码为:

1
2
3
4
5
Ff = fft(f_);
Ff(abs(Ff)<5)=0.0; %将绝对值小于5的傅里叶系数都设为0.0
f1 = ifft(Ff);
plot(x, f_, x, f1);
legend('带噪声函数', '离散傅里叶变换结果');

离散傅里叶变换去噪的局限性

然而,如果将三角函数进行阶梯函数叠加使其成为分片光滑函数,傅里叶变换结合硬阈值(去掉绝对值小于𝝀的系数)算法的效果,要么无法消除噪声,要么会破坏信号信息,导致恢复效果甚至明显不如五点均值的效果

软阈值算和硬阈值的区别

  • 硬阈值算子 $T_{\lambda}(y)=\left\{\begin{array}{l}y,|y| \geq \lambda \\ 0,|y|<\lambda\end{array},\right.$ 注意到其绝对值不小于$\lambda$的元素均末发生变化。对y的傅里早系数做硬阈值算子处理相当于求解问题$\operatorname{min}$ $_{x} \frac{1}{2}|x-y|_{2}^{2}+\frac{1}{2} \lambda^{2}|\widehat{x}|_{0},$ 其中, $\widehat{x}$ 表示 $x$ 的离散傅里叶变换, $|\boldsymbol{t}|_{0}$ 表示 $\boldsymbol{t}$ 的非零元素的个数。

  • 软阈值算子 $T_{\lambda}(y)=\operatorname{sign}(y) \cdot \max (|y|-\lambda, 0)$ ,对 $y$ 的傅里叶系数做软阈值算子处理相当于求解问题$\min _{x} \frac{1}{2}|x-y|_{2}^{2}+\lambda|\hat{x}|_{1}$

软阈值算法完成高斯噪声去噪

对于前例,软阈值算法可以获得比硬阈值算法更好的效果,但是仍然不如五点均值滤波,原因是函数并不具备全局光滑性或周期性,导致傅里叶系数的稀疏度不足。代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
clear,close all,rng default
x=0:pi/50:2*pi;f = sin(x);
for i = 1:5 %强行设置若干处阶梯越阶,构造分片光滑函数
f(i*20-19:i*20)=f(i*20-19:i*20)+(5-i);
end
y = f + randn(size(x))*0.5;
Ff = fft(y);
Ff = sign(Ff).*max(abs(Ff)-5,0);
%将绝对值小于5的傅里叶系数都设为0.0,模大于5的系数幅角不变,模减5

f1 = ifft(Ff);
snr(f,f1-f) %结果的SNR=16.9631

总结

数字信号处理的去噪问题往往种类繁多,最佳的去噪方法也往往不是唯一的,既与噪声类型和强度有关,也与目标函数(真实信号值)的函数特性有关。

  • 常见的噪声类型除高斯噪声外,还有冲击噪声(随机的部分采样值出现与真实信号值无关的错误函数值),通常可以使用中值滤波器进行噪点判别或直接调整。

  • 泊松噪声是信号处理较难解决的非加性噪声之一,在一个较大的假想最大值的基础上,在每点设$\lambda = k \cdot$ 真实值,$k>0$为比例因子,然后基于参数(期望)为$\lambda$的泊松分布下获取一个整数值$Y$,带噪信号的该点对应信号值为$\frac Yk$。容易发现$k$越小时,泊松噪声造成的影响越大。

  • 傅里叶变换对于三角函数去噪效果极佳,因为三角函数在此变换系数满足完美稀疏,对于分片光滑函数效果则一般,根据目标函数特性选择合适变换及约束(正则化)方法很关键

声音信号处理

  • 读取音乐的函数为audioinfo,返回值为yFs,其中Fs为一秒钟内信号读取次数,在这里(大多数音乐文件中)Fs为14400
  • 播放音乐的函数为soundsound函数中可带参数Fs,可以利用这个参数调节音乐的频率,如下例中将Fs除以$\sqrt{2}$后会发现音乐明显频率降低了。
1
2
3
4
5
6
info = audioinfo('Week14音乐识别文件库\Music1.mp3');
[y, Fs]= audioread('Week14音乐识别文件库\Music1.mp3');
ys = y(1:11*Fs, 1);
sound(ys, Fs);
sound(ys, Fs/sqrt(2));
plot(abs(fft(ys)));
  • 版权声明: 本博客所有文章除特别声明外,著作权归作者所有。转载请注明出处!

扫一扫,分享到微信

微信分享二维码
  • Copyrights © 2020-2021 chenk
  • 由 帅气的CK本尊 强力驱动
  • 访问人数: | 浏览次数:

请我喝杯咖啡吧~

支付宝
微信