搜索
１／ｆ噪声的仿真产生和验证SimulationandVerificationof1/fNoiseWUYuan-zhen，WUMU-zhiSuzhouTestingInstrumentFactorySuzhouJiangsu215129China:Inthispaper,theprinciplesofthe1/fnoisearediscussed.Inaccordancewiththeseprinciples,usingthewhitenoisetobeaddedtosimplifytheright1/fdigitalfiltermethod;theauthorwroteprogramsbasedonMatlabandgetsabettersimulationof1/fnoise.Andthroughthecommoncomputingpowerspectraldensitywayofthealgorithm,theverificationisgeneratedbythe1/fnoiseperformance,the1/fnoisecharacteristicsareshown.近年来，随着在物理、化学、生物医学等基础学科中对自然现象及规律的深入探索，在通信、测量、导航、医疗仪器等工业领域中人们必须测量及接受越来越小的电信号。测量或者接受装置中的内部噪声影响越来越为明显，这种噪声已经成为制约电子仪器和装备进一步提高信号检测灵敏度及改进接受信号质量的关键因素。仿真产生噪声信号可以帮助研究人员研究克服噪声的方法，帮助技术人员研制出高抗干扰的仪器设备或指出改进产品的方向。1/f噪声广泛存在于电子元器件和电子线路中，从1/f噪声测量分析中可以得到很多有用的信息。本文主要讨论1/f噪声的特点和仿真产生的方法。Matlab的信号处理工具箱是信号算法文件的集合，它处理的对象是信号和系统，信号处理工具箱提供了很多命令行函数来进行谱分析，包括经典的无界技术以及特征值技术。另外，也添加了很多对象，用以提高可用性。因此我们使用该工具箱来仿真产生1/f噪声源，并且用其来分析验证仿真产生出的噪声可信度。11/f噪声的特性1.11/f噪声（低频率噪声）的特性电子器件中的1/f噪声具有两个基本特性：1)1/f噪声在一个相当宽的频带内，功率谱密度与频率f成反比,且频带上、下限都是有限值，上限频率视1/f噪声与背景白噪声的相对大小而定；下限频率接近直流时功率密度仍呈很好的1/f特性。2)1/f噪声电压和电流的功率谱密度近似与流经器件的电流的平方成正比，这意味着1/f噪声起源于电阻的涨落。设流过电阻R的电流保持恒定，电阻存在着涨落，引起电阻两端电压涨落，则有:δV=1*δR(1-1)Sv(f)=I2SR(f)I2(1-2)假设功率谱密度是Pnf=Sf/f，相应的电压密度是e2nf=Af/f。从频率f1到f2之间的总噪声是v2nf=f2f1Affdf=Aflog(f2f1)(1-3)所以，1/f谱图取决于上下的截止频率，而不是绝对带宽。1.21/f噪声的自相关函数LTV(线性时不变)系统的输出可以用卷积积分来表示：(1-4)式中t是输出的观察时间。通常我们更关心的是计算系统输出的均方根值，假设输入是由一个静态噪声源驱动的，为了得到这个结果，第一步是计算输出均方根值x~2=?馈?-∞?馈?-∞h(u)h(v)m(t-u)m(t-v)ni(t-u)ni(t-v）dudv(1-5)这里f(t)被一个噪声信号n(t)取代。均方值通过计算整个时间段的平均值可以得到：x-2=limk→∞1k?馈?-∞?馈?-∞h(u)h(v)m(t-u)m(t-v)ni(t-u)ni(t-v)dudvx-2=?馈?-∞?馈?-∞h(u)h(v)m(t-u)m(t-v)［limk→∞1kni(t-u)ni(t-v)］dudv(1-7)ni(t)表示的是i次抽样函数。式(1-7)中间括号的内容是随机输入噪声信号的自相关函数，所以可以简化成为x-2?馈?-∞?馈?-∞h(u)h(v)m(t-u)m(t-v)［Rnn(u-v)］dudv（1-8）上面这个式子指出了自相关函数在时域噪声分析中的重要性。在噪声处理中最常见的两种噪声是白噪声和1/f噪声的自相关函数。理想的1/f噪声是一个多样的过程，它的自相关函数无法用闭合的形式来简单的表达出来。实际应用中进行了一些近似。白噪声的双边功率谱密度(double-sidedPSD)是v2Δf=N02（1-9）N0是单边噪声功率谱密度，单位是V2/HZ。自相关函数就是简单的功率谱密度的富利叶反变换。所以，白噪声的自相关函数由冲击响应函数给出Rnn(τ)=N02δ(τ)(1-10)1/f噪声，可以看成是白噪声通过一个噪声滤波器得到的结果。白噪声驱动的通过线性时间滤波器的输出的自相关函数是Ryy(τ)=N02?馈?-∞h(λ)h(τ+λ)dλ(1-11)h(t)是成型滤波器的冲击响应。理想的1/f成型噪声滤波器得冲击响应是h(t)2fct,t＞0(1-11)fc是转角频率。1/f噪声的自相关函数从下式计算得到：Ryy(τ)=N0fc?馈蕺?01λ1λ+τdλ(1-13)可以得到Ryy=2N0fc1n［λ+λ+τ］∞0(1-14)上面这个式子形式是非常简单的，可惜的是，式(1-14)是发散的，所以并不可能去精确的计算得到1/f噪声的自相关函数。我们可以简化得到有理形式的滤波器（或者一串滤波器），解决方法是定义一个总的操作时间，并使用恰当的近似，式(1-14)可以转变为Ryy(top,τ)=2N0fc1n［λ+λ+τ］top-τδ(1-15)假设top?恙营恙模?可以得到Ryy(top,τ)≌N0fc［1n(4)+1n(topfc）-1n(τfc)］（1-16）注意，top和δ，是等同于时域内的两个频率限fmin(等同于/top)和(fc等同于1/δ)，所以，整个1/f噪声过程包括了N个10倍频段，(1-12)可以写成N=log10(topfc)（1-17）Ryy(N,τ)≌N0fc［1n(4)+2.3N-1n(τfc)］(1-18)在1/f噪声仿真的实际应用中，式(1-14)比式(1-19)更实用。1.31/f噪声的产生方法产生1/f噪声的最直观的方法就是将白噪声通过噪声滤波器过滤之后得到。但是，理想的1/f滤波器并不是有理的，所以需要开发一个可以近似得到1/f噪声的滤波器。需解决的主要问题是求这一近似的滤波器的算法。使用下图所示简化的噪声滤波器。图1仿真1/f噪声成型滤波器的示意图一系列的hn(t)进行叠加，可以得到一系列的滤波器，产生了完整的1/f噪声的输出，传输函数如下给出H(s)=1+∑N+1n=010-0.5n22πfcs+2π10-nfc(1-19)N定义了总频率带宽的10倍频数量，fc是转角频率。1/f噪声的数字化产生需要将模拟滤波器变成数字化滤波器，z域传输函数很容易计算得到，计算的结果如下式：H(z)=1+∑Nn=010-0.5n2πfcTs1-z-1exp(-2π10-nfcTs)(1-20)现在可以使用Matlab内置的filter函数来进行处理了。在仿真产生1/f噪声时，比较合理N的值大约为10。这意味着一个带宽为1MHz的系统，最低频率在10-4左右。这样，总的仿真时间长度需要在s。一般来说104s，仿真N个10倍频段需要大约10N个点。如果N6，考虑到对计算机内存和资源的要求，是不现实的。在实际应用中，通常选定采样频率为107，模拟时间长度为10-3。考虑到所需要的精度大概需要104-105个点。为了解决这个问题我们定义一个参量NeffNeff=1+log10(tsimTs)(1-21)tsim是总的仿真时长。然后整个数字滤波器被分为两个部分，如图2所示，所有的nNeff子滤波器，按照式(1-20)来计算，所有的n＞Neff部分，输出在整个仿真过程中几乎都是常量。意味着高阶的子滤波器可以进行简单的将数值（偏移量）相加。每一部分所得均方根等于稳态时子滤波器的输出均方根。图2数字滤波器可以被Neff分为两个部分1.4产生白噪声并计算其功率谱密度1.4.1使用matlab仿真产生白噪声如上节分析，因为仿真产生1/f噪声的方法是白噪声通过一个滤波器，所以先需要仿真得到一个白噪声。(程序较简单略)在Matlab中产生白噪声，可以使用randn函数的得到一个具有高斯正态分布的白噪声noise=rms×randn(1,npts)产生得到一个均值为0，均方根值为rms，总点数为npts的噪声，均方根值由白噪声底阶和采样频率共同决定，如下式得到：rms=N02Ts(1-22)N0是白噪声的单边功率谱密度，Ts是采样频率。采样频率最好取为10的因子并且略大于系统带宽，这样在整个带宽内白噪声都有一个平坦的功率谱。程序输入输出量说明：functionnoise=shot_noise(Ts,N0,npts)Ts采样时间，fs=1/Ts采样频率，N0白噪声的单边功率谱密度，npts仿真所取的点，noise为所得到的噪声数值序列图3仿真得到的白噪声时域波形1.4.2计算功率谱密度的三种不同的方法1.4.3白噪声的功率谱密度为了验证产生的白噪声的性能，计算它的PSD进行验证functionpsd_s(Ts,N0,npts,method)%PSDofshotnoise除了method外与产生程序的变量意义相同，method为产生功率谱密度算法参数。图4仿真得到的白噪声功率谱1.5仿真产生1/f噪声functionnoise=f1_noise(tsim,Ts,fc,N0,N,npts)%noise=f1_noise(tsim,Ts,fc,N0,N,npts)%variabledeclaration%Ts:采样时间长度%N:10倍频数辆%tsim:总仿真时长%fs:采样周期%fc:转角频率具体源程序代码见2.1.1节图5仿真得到的1/f噪声时域波形tsim,仿真持续时间，Ts采样时间，N10倍频段，fs采样频率，npts总共仿真的点数。得到了一个杂乱无章的1/f噪声时域波型，需要验证它就是我们所需要的。1.61/f噪声的功率谱密度functionpsd_f(tsim,Ts,fc,N0,N,npts,method)%pad_f(tsim,Ts,fc,N0,N,npts,method)%variabledeclaration%Ts:SamplingPeriod%N:frequecncydecades%tsim:totalsimulatedtimeinterval%fs:SamplingFrequency%fc:cornerfrequncy具体的源程序代码见2.1.2节。除了method外与产生程序的变量意义相同，method为产生功率谱密度算法参数。direct用直接法进行计算，直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。ndirect，间接法由序列x(n)估计出自相关函数R(n)，对R(n)进行傅立叶变换，得到x(n)的功率谱估计。其中complex方法，会将加入三种窗函数的图谱估计直接显示在一幅图上。图6仿真得到的1/f噪声功率谱密度可以看到，PSD近似与lnf成一个反比的关系，随着f的增大而减小。对图形进行处理，横轴以ln（1/f）为坐标，下降的关系是成直线的，在低频率范围内出现的扭曲，是滤波器传递函数在N＞Neff时的近似造成的。图7对数坐标-1/f功率谱密度取对数坐标作图(见图7)，舍弃超过转角频率的点，进行线性拟合得到：b=-1.4486a=-2.2333r=-3.99222仿真过程的实现2.1仿真产生1/f噪声2.1.1仿真产生1/f噪声的程序functionnoise=f1_noise(tsim,Ts,fc,N0,N,npts)%noise=f1_noise(tsim,Ts,fc,N0,N,npts)%v