本文资料及代码来源:http://wiki.woodpecker.org.cn/moin/MiscItems/2011-08-24
以下是自己的相关整理,代码及文章内,请支持上述原创!
- 引言: matlab中有agwn(signal,SNR)函数可以再原始信号的基础上增加SNR信噪比的噪音;但貌似python的若干库中都没有现成的代码?那么怎么办呢?
- 答案:理解信噪比的定义,手动写一个呗!
定义
原始信号:x
噪声信号:n
信噪比:SNR
SNR=10log10PsignalPnoise=10log10∑x2∑n2
wgn(x, snr)
中x为信号,snr为信噪比,返回满足条件的高斯白噪声,只需要:
x += wgn(x, snr)
即可以得到和matlab的awgn相同的效果。
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Jun 23 18:08:20 2017
@author: brucelau
"""
import numpy as np
import matplotlib.pyplot as plt
r = np.linspace(0.000001,20)
snr = 10*np.log10(r)
plt.plot(r,snr)
plt.grid()
#plt.xlabel('$\farc{{\sum{x^2}}{\sum{n^2}}}$')
plt.xlabel('$10*log_{10}\\frac{\sum{x^2}}{\sum{n^2}}$')
plt.ylabel('SNR')
plt.title('The ratio-snr curve')
plt.tight_layout()
结果
代码
# -*- coding: utf-8 -*-
"""
Created on Mon Jun 19 21:21:58 2017
@author: Administrator
"""
'''
程序中用hist()检查噪声是否是高斯分布,psd()检查功率谱密度是否为常数。
'''
import numpy as np
import pylab as plt
def wgn(x, snr):
snr = 10**(snr/10.0)
xpower = np.sum(x**2)/len(x)
npower = xpower / snr
return np.random.randn(len(x)) * np.sqrt(npower)
plt.figure()
t = np.arange(0, 10000000) * 0.1
x = np.sin(t)
n = wgn(x, 6)
xn = x+n # 增加了6dBz信噪比噪声的信号
plt.subplot(221)
plt.plot(t[:1024],x[:1024])
plt.title('The original signal-x')
plt.subplot(222)
plt.plot(t[:1024],xn[:1024])
plt.title('The original sinal with Gauss White Noise')
plt.subplot(223)
plt.hist(n, bins=100, normed=True)
plt.title('Gauss Noise Distribution')
plt.subplot(224)
plt.psd(n)
plt.title('PSD')
plt.tight_layout()
plt.savefig('show the result.png',dpi=600)