数值微分的精度

numpy_diff.py

本习题的程序文件

在系统控制和信号分析中,经常会需要对信号进行微分运算。根据微分的数学定义,只需要将信号的连续两个数值的差除以取样时间间隔既可计算出信号的微分。然而实际系统中的随机噪声和量化噪声对微分的结果有显著的影响,本习题对数值微分的精度问题进行探讨,找到特定量化噪声情况下的最佳数值微分方案。

量化噪声

将模拟信号转换成数字信号时,会产生量化噪声。这是因为我们无法用有限的比特数表示无限多种可能的模拟信号。例如若用一个16比特的整数表示0V到10V之间变化的电压信号,那么数字信号所能表示的最小的电压间隔为:10/2^{16},约等于0.00015V。

两点微分

函数f的微分可以用下面的公式计算,只要取足够小的h值,就可以使用f(x)f(x-h)计算f'(x)

f'(x) = \lim_{h\to0} \frac{f(x) - f(x-h)}{h}

请根据上面的微分公式编写程序计算数组y的微分:

import numpy as np

def err(x, y):
    return np.average(abs(x-y))

h = 0.0001
x = np.arange(0, 2*np.pi, h)
f = np.sin(x)  # 需要进行微分的数据
df = np.cos(x) # 理想的微分值

【你的程序】用二点公式计算f的微分df2

print "two point differentiation"
print err(df[:-1], df2)
print err(df[1:],  df2)
用numpy.diff(a)可以计算数组a中前后两个数值的差。

微分之后的数组长度为原数组的长度减1,因此使用“df[:-1]”或“df[1:]”与数值微分的结果计算误差。程序的输出应为:

two point differentiation
3.18314207566e-05
3.18314207755e-05

由结果可知使用f(x+h)f(x)计算的微分与使用f(x)f(x-h)计算的微分的误差是相同的。

多点微分

进行数值微分时,使用的点数越多,则微分值越精确。可以从下面的链接查找到不同点数的微分系数。

本习题采用此微分系数表中的后向微分系数计算,即采用f(x)以及f(x-k \cdot h)计算f'(x)。在数据采集系统中,只有后向微分能实时地进行计算,因为时刻t的微分只需要用t之前的数据计算。例如,使用f(x), f(x-h), f(x-2h)计算f'(x)的公式如下:

f'(x) = \frac{3f(x) -4f(x-h) + f(x-2h)}{2h} + O(h^2)

其中O(h^2)为误差项,它和h^2成正比。由于h是一个非常小的值,O(h^2)比两点微分的误差项O(h)要小很多。下面请编写程序计算三点微分:

【你的程序】用三点公式计算f的微分df3

print "three point differentiation"
print err(df[2:], df3)

程序的输出为:

three point differentiation
2.12203334556e-09

可以看出对于同样的h,三点微分的精度要比二点微分高许多。这意味着我们可以适当选择较大的h进行运算,而用较大的h计算微分能够减低噪声的影响。

最后,请编写能计算N点微分的函数diff_npoint():

from numpy.lib.stride_tricks import as_strided

def diff_npoint(a, coeffs, h):
    "计算数组a的任意点的微分,coeffs为f(x-k*h),...,f(x-h),f(x)的系数"
    【你的程序】

print "four point differentiation"
print err(df[3:], diff_npoint(f, [-1.0/3, 1.5, -3, 11.0/6], h))

程序的输出应为:

four point differentiation
4.07252437767e-12
在diff_npoint()之内可以使用循环做求和运算,但为了显得更酷一些,你能用as_strided()写出无循环的程序吗?

减小量化噪声的影响

如果要微分的信号是用模数转换器采样说得的电压信号,那么由于模数转换器的精度,信号中存在量化噪声,请编写程序将数组f转换为1024量化等级的数组fq。fq的最小值和最大值仍然为-1和1,但是它只包含1024个不同的数值:

【你的程序】将f量化为4096个等级的数组fq

print "three point differentiation of quantized signal"
print err(df[2:], diff_npoint(fq, [0.5, -2, 1.5], h))
import pylab as pl
pl.plot(x,f,x,fq,lw=2)
/tech/static/books/scipyex/_images//numpy_diff_01.png

模拟信号和量化之后的数字信号

【图:模拟信号和量化之后的数字信号】为数组f和fq的局部放大,可以看出fq呈阶梯状。直接计算fq的微分会出现很大的量化噪声。程序的输出为:

three point differentiation of quantized signal
1.80630449261

为了减小量化噪声,可以将公式中的h替换为k*h,增加前后两个取样值的X轴的距离:

f'(x) = \frac{3f(x) -4f(x-k\cdot h) + f(x-2k\cdot h)}{2k\cdot h}

请编写diffq_3point(),用带k的三点微分公式计算数值微分:

def diffq_3point(a, k, h):
    "用3点微分公式计算数组a的微分"
    【你的程序】

ks = range(20,2000,20) #
errors = [err(df[2*k:], diffq_3point(fq,k,h)) for k in ks] #
pl.figure()
pl.plot(ks, errors, lw=2)
pl.xlabel("k")
pl.ylabel("error")

❶让k值在20到2000之间变化,❷并计算每个k值所对应的数值微分和理论值之间的误差。【图:k值和误差之间的关系】显示了误差和k之间的关系,由图可以看出对于量化信号存在一个最优的k值使得其微分的误差最小。

/tech/static/books/scipyex/_images//numpy_diff_02.png

k值和误差之间的关系

请读者试验当数据中存在随机噪声时,k值和误差之间的关系。

內容目录

上一个主题

分析电机编码器的A相和Z相波形

下一个主题

拟合正弦波的频率、相位和振幅

本页

loading...