摘 要：本文介绍了1/f噪声的含义和原理，依照这些原理，通过白噪声加简化得到的有理1/f数字滤波器的方法，使用Matlab编程，仿真产生了真实1/f噪声。并根据常见的计算功率谱密度的算法，验证了仿真产生的1/f噪声的性能。
　　关键词:噪声；1/f噪声；仿真；功率谱密度；
　　中图法分类号:TB533 文献标志码:B
　　
　　Simulation and Verification of 1/f Noise
　　WU Yuan-zhen，WU MU-zhi
　　Suzhou Testing Instrument Factory Suzhou Jiangsu 215129 China
　　ABSTRACT:In this paper, the principles of the 1/f noise are discussed. In accordance with these principles, using the white noise to be added to simplify the right 1/f digital filter method; the author wrote programs based on Matlab and gets a better simulation of 1/f noise. And through the common computing power spectral density way of the algorithm, the verification is generated by the 1/f noise performance, the 1/f noise characteristics are shown.
　　Key words:Matlab;noise;1/f noise;simulation;PSD
　　
　　近年来，随着在物理、化学、生物医学等基础学科中对自然现象及规律的深入探索，在通信、测量、导航、医疗仪器等工业领域中人们必须测量及接受越来越小的电信号。测量或者接受装置中的内部噪声影响越来越为明显，这种噪声已经成为制约电子仪器和装备进一步提高信号检测灵敏度及改进接受信号质量的关键因素。仿真产生噪声信号可以帮助研究人员研究克服噪声的方法，帮助技术人员研制出高抗干扰的仪器设备或指出改进产品的方向。
　　1/f噪声广泛存在于电子元器件和电子线路中，从1/f噪声测量分析中可以得到很多有用的信息。本文主要讨论1/f噪声的特点和仿真产生的方法。
　　Matlab的信号处理工具箱是信号算法文件的集合，它处理的对象是信号和系统，信号处理工具箱提供了很多命令行函数来进行谱分析，包括经典的无界技术以及特征值技术。另外，也添加了很多对象，用以提高可用性。因此我们使用该工具箱来仿真产生1/f噪声源，并且用其来分析验证仿真产生出的噪声可信度。
　　1 1/f噪声的特性
　　1.1 1/f噪声（低频率噪声）的特性
　　电子器件中的1/f 噪声具有两个基本特性：
　　1)1/f 噪声在一个相当宽的频带内，功率谱密度与频率f 成反比, 且频带上、下限都是有限值，上限频率视1/f 噪声与背景白噪声的相对大小而定；下限频率接近直流时功率密度仍呈很好的1/f 特性。
　　2)1/f 噪声电压和电流的功率谱密度近似与流经器件的电流的平方成正比，这意味着1/f噪声起源于电阻的涨落。设流过电阻R的电流保持恒定，电阻存在着涨落，引起电阻两端电压涨落，则有:
　　 δV=1*δR(1-1)
　　 Sv(f)=I2SR(f)I2(1-2)
　　假设功率谱密度是Pnf=Sf/f，相应的电压密度是e2nf=Af/f。从频率f1到f2之间的总噪声是
　　v2nf=f2f1Affdf=Aflog(f2f1)(1-3)
　　所以，1/f谱图取决于上下的截止频率，而不是绝对带宽。
　　1.2 1/f噪声的自相关函数
　　LTV(线性时不变)系统的输出可以用卷积积分来表示：
　　(1-4)
　　式中t是输出的观察时间。
　　通常我们更关心的是计算系统输出的均方根值，假设输入是由一个静态噪声源驱动的，为了得到这个结果，第一步是计算输出均方根值
　　x~2=∞-∞∞-∞h(u)h(v)m(t-u)m(t-v)ni(t-u)ni(t-v）dudv(1-5)
　　这里f(t)被一个噪声信号n(t)取代。均方值通过计算整个时间段的平均值可以得到：
　　x-2=limk→∞1k∞-∞∞-∞h(u)h(v)m(t-u)m(t-v)ni(t-u)ni(t-v)dudv
　　x-2=∞-∞∞-∞h(u)h(v)m(t-u)m(t-v)［limk→∞1kni(t-u)ni(t-v)］dudv(1-7)
　　ni(t)表示的是i次抽样函数。
　　式(1-7)中间括号的内容是随机输入噪声信号的自相关函数，所以可以简化成为
　　x-2∞-∞∞-∞h(u)h(v)m(t-u)m(t-v)［Rnn(u-v)］dudv（1-8）
　　上面这个式子指出了自相关函数在时域噪声分析中的重要性。在噪声处理中最常见的两种噪声是白噪声和1/f噪声的自相关函数。理想的1/f噪声是一个多样的过程，它的自相关函数无法用闭合的形式来简单的表达出来。实际应用中进行了一些近似。
　　白噪声的双边功率谱密度(double-sidedPSD)是
　　 v2Δf=N02（1-9）
　　N0是单边噪声功率谱密度，单位是V2/HZ。自相关函数就是简单的功率谱密度的富利叶反变换。所以，白噪声的自相关函数由冲击响应函数给出
　　 Rnn(τ)=N02δ(τ)(1-10)
　　1/f噪声，可以看成是白噪声通过一个噪声滤波器得到的结果。白噪声驱动的通过线性时间滤波器的输出的自相关函数是
　　Ryy(τ)=N02∞-∞h(λ)h(τ+λ)dλ(1-11)
　　h(t)是成型滤波器的冲击响应。理想的1/f成型噪声滤波器得冲击响应是
　　h(t)2fct,t＞0(1-11)
　　fc是转角频率。1/f噪声的自相关函数从下式计算得到：
　　Ryy(τ)=N0fc∞01λ1λ+τdλ(1-13)
　　可以得到
　　Ryy=2N0fc1n［λ+λ+τ］∞0(1-14)
　　上面这个式子形式是非常简单的，可惜的是，式(1-14)是发散的，所以并不可能去精确的计算得到1/f噪声的自相关函数。我们可以简化得到有理形式的滤波器（或者一串滤波器），解决方法是定义一个总的操作时间，并使用恰当的近似，式(1-14)可以转变为
　　Ryy(top,τ)=2N0fc1n［λ+λ+τ］top-τδ(1-15)
　　假设topτδ，可以得到
　　Ryy(top,τ)≌N0fc［1n(4)+1n(topfc）-1n(τfc) ］ （1-16）
　　注意，top和δ，是等同于时域内的两个频率限fmin(等同于/top)和(fc等同于1/δ)，所以，整个1/f噪声过程包括了N个10倍频段，(1-12)可以写成
　　N=log10(topfc)（1-17）
　　Ryy(N,τ)≌N0fc［1n(4)+2.3N-1n(τfc)］(1-18)
　　在1/f噪声仿真的实际应用中，式(1-14)比式(1-19)更实用。
　　1.3 1/f噪声的产生方法
　　产生1/f噪声的最直观的方法就是将白噪声通过噪声滤波器过滤之后得到。但是，理想的1/f滤波器并不是有理的，所以需要开发一个可以近似得到1/f噪声的滤波器。需解决的主要问题是求这一近似的滤波器的算法。
　　使用下图所示简化的噪声滤波器。
　　图1 仿真1/f噪声成型滤波器的示意图
　　一系列的hn(t)进行叠加，可以得到一系列的滤波器，产生了完整的1/f噪声的输出，传输函数如下给出
　　H(s)=1+∑N+1n=010-0.5n22πfcs+2π10-nfc(1-19)
　　N定义了总频率带宽的10倍频数量，fc是转角频率。1/f噪声的数字化产生需要将模拟滤波器变成数字化滤波器，z域传输函数很容易计算得到，计算的结果如下式：
　　H(z)=1+∑Nn=010-0.5n2πfcTs1-z-1exp(-2π10-nfcTs)(1-20)
　　现在可以使用Matlab内置的filter函数来进行处理了。在仿真产生1/f噪声时，比较合理N的值大约为10。这意味着一个带宽为1MHz的系统，最低频率在10-4左右。这样，总的仿真时间长度需要在s。一般来说104s，仿真N个10倍频段需要大约10N个点。如果N>6，考虑到对计算机内存和资源的要求，是不现实的。在实际应用中，通常选定采样频率为107，模拟时间长度为10-3。考虑到所需要的精度大概需要104-105个点。为了解决这个问题我们定义一个参量Neff
　　Neff=1+log10(tsimTs)(1-21)
　　tsim是总的仿真时长。然后整个数字滤波器被分为两个部分，如图2所示，所有的nNeff子滤波器，按照式(1-20)来计算，所有的n＞Neff部分，输出在整个仿真过程中几乎都是常量。意味着高阶的子滤波器可以进行简单的将数值（偏移量）相加。每一部分所得均方根等于稳态时子滤波器的输出均方根。
　　图2 数字滤波器可以被Neff分为两个部分
　　1.4 产生白噪声并计算其功率谱密度
　　1.4.1 使用matlab仿真产生白噪声
　　如上节分析，因为仿真产生1/f噪声的方法是白噪声通过一个滤波器，所以先需要仿真得到一个白噪声。(程序较简单略)
　　在Matlab中产生白噪声，可以使用randn函数的得到一个具有高斯正态分布的白噪声
　　noise=rms×randn(1,npts)
　　产生得到一个均值为0，均方根值为rms，总点数为npts的噪声，均方根值由白噪声底阶和采样频率共同决定，如下式得到：
　　 rms=N02Ts(1-22)
　　N0是白噪声的单边功率谱密度，Ts是采样频率。采样频率最好取为10的因子并且略大于系统带宽，这样在整个带宽内白噪声都有一个平坦的功率谱。
　　程序输入输出量说明：
　　function noise=shot_noise(Ts,N0,npts)
　　Ts采样时间，fs=1/Ts采样频率，N0白噪声的单边功率谱密度，npts仿真所取的点，noise为所得到的噪声数值序列
　　图3 仿真得到的白噪声时域波形
　　1.4.2 计算功率谱密度的三种不同的方法
　　函数中使用direct关键词则使用了直接法进行计算，又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。
　　使用indirect关键词，使用间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，得到x(n)的功率谱估计。
　　其中complex关键词，综合方法会将加入三种窗函数(Hamming,Kaiser,Chebyshev)的图谱估计直接显示在一幅图上，比较直观和明显。
　　1.4.3 白噪声的功率谱密度
　　为了验证产生的白噪声的性能，计算它的PSD进行验证
　　function psd_s(Ts,N0,npts,method)
　　 % PSD of shot noise
　　除了method外与产生程序的变量意义相同，method为产生功率谱密度算法参数。
　　图4 仿真得到的白噪声功率谱
　　1.5 仿真产生1/f噪声
　　function noise=f1_noise(tsim,Ts,fc,N0,N,npts)
　　% noise=f1_noise(tsim,Ts,fc,N0,N,npts)
　　% variable declaration
　　% Ts: 采样时间长度
　　% N:10倍频数辆
　　% tsim:总仿真时长
　　% fs:采样周期
　　% fc:转角频率
　　具体源程序代码见2.1.1节
　　图5 仿真得到的1/f噪声时域波形
　　tsim,仿真持续时间，Ts采样时间，N10倍频段，fs采样频率，npts总共仿真的点数。
　　得到了一个杂乱无章的1/f噪声时域波型，需要验证它就是我们所需要的。
　　1.6 1/f噪声的功率谱密度
　　function psd_f(tsim,Ts,fc,N0,N,npts,method)
　　% pad_f(tsim,Ts,fc,N0,N,npts,method)
　　% variable declaration
　　% Ts: Sampling Period
　　% N:frequecncy decades
　　% tsim:total simulated time interval
　　% fs:Sampling Frequency
　　% fc:corner frequncy
　　具体的源程序代码见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.4486 a=-2.2333r=-3.9922
　　2 仿真过程的实现
　　2.1 仿真产生1/f噪声
　　2.1.1 仿真产生1/f噪声的程序
　　function noise=f1_noise(tsim,Ts,fc,N0,N,npts)
　　% noise=f1_noise(tsim,Ts,fc,N0,N,npts)
　　% variable declaration
　　% Ts: Sampling Period
　　% N:frequecncy decades
　　% tsim:total simulated time interval
　　% fs:Sampling Frequency
　　% fc:corner frequncy
　　Neff=1+log10(tsim/Ts);%parameter Neff
　　noise_tmp=shot_noise(Ts,N0,npts);
　　fs=1/Ts;
　　% Generation white noise using the shot_noise function
　　fn=0;
　　for n=0:N
　　if n<=Neff
　　b=［0 ((10^(-0.5*n))*sqrt(2)*pi*fc/fs)］;
　　a=［1-1*exp(-2*pi*(10^(-n))*fc/fs)］;
　　fn=fn+filter(b,a,noise_tmp);
　　else
　　fn=fn+sqrt(N0*pi*fc/4)*randn(1,1);
　　 end
　　end
　　plot(fn);pause;
　　noise=fn;
　　2.1.2 1/f噪声功率谱密度计算程序
　　function psd_f(tsim,Ts,fc,N0,N,npts,method)
　　% pad_f(tsim,Ts,fc,N0,N,npts,method)
　　% variable declaration
　　%Ts: Sampling Period
　　%N:frequecncy decades
　　%tsim: total simulated time interval
　　%fs: Sampling Frequency
　　%fc: corner frequncy
　　% methods
　　% 'direct'
　　% 'indirect'
　　% 'periodogram'
　　% 'Kaiser'
　　% 'Chebyshev'
　　% 'complex'
　　noise=f1_noise(tsim,Ts,fc,N0,N,npts);
　　Fs=1/Ts;
　　n=N0;
　　xn=noise;
　　%direct method
　　switch lower(method)
　　 case{'direct'}
　　 window= boxcar(length(noise));
　　 nfft=1024;
　　 ［Pxx,f］=periodogram(noise,window,nfft,Fs);
　　 plot(f,10*log10(Pxx))
　　%indirect method
　　case{'indirect'}
　　 nfft=1024;
　　 window=boxcar(length(n));
　　 noverlap=0;
　　 p=0.9;
　　 ［Pxx,Pxxc］=psd(xn,nfft,Fs,window,noverlap,p);
　　 index=0:round(nfft/2-1);
　　 k=index*Fs/nfft;
　　plot_Pxx=10*log10(Pxx(index+1));
　　 plot_Pxxc=10*log10(Pxxc(index+1));
　　 figure(1)
　　 plot(k,plot_Pxx);
　　 pause;
　　 figure(2)
　　plot(k,［plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc］);
　　case {'periodogram'}
　　 % Instantiate spectrum object and call its PSD method.
　　h=spectrum.periodogram('rectangular');
　　hopts=psdopts(h,xn);% Default options based on the signal x
　　set(hopts,'Fs',Fs,'SpectrumType','twosided','CenterDC',true);
　　 psd(h,xn,hopts)
　　 case{'kaiser'}
　　 h=spectrum.welch('hamming',64);
　　hpsd1=psd(h,xn,'Fs',Fs);
　　W=hpsd1.Frequencies;
　　h.WindowName='Kaiser';
　　hpsd2=psd(h,xn,'Fs',Fs);
　　Pxx2=hpsd2.Data;
　　hpsd=dspdata.psd(Pxx2,W,'Fs',Fs);
　　plot(hpsd);
　　legend('kaiser');
　　case{'chebyshev'}
　　h=spectrum.welch('hamming',64);
　　hpsd1=psd(h,xn,'Fs',Fs);
　　W=hpsd1.Frequencies;
　　h.WindowName='Chebyshev';
　　hpsd3=psd(h,xn,'Fs',Fs);
　　Pxx3=hpsd3.Data;
　　hpsd=dspdata.psd(Pxx3,W,'Fs',Fs);
　　plot(hpsd);
　　legend('Chebyshev');
　　case{'complex'}
　　h= spectrum.welch('hamming',64);
　　% Create three power spectral density estimates.
　　hpsd1=psd(h,xn,'Fs',Fs);
　　Pxx1=hpsd1.Data;
　　W=hpsd1.Frequencies;
　　h.WindowName='Kaiser';
　　hpsd2=psd(h,xn,'Fs',Fs);
　　Pxx2=hpsd2.Data;
　　h.WindowName='Chebyshev';
　　hpsd3=psd(h,xn,'Fs',Fs);
　　Pxx3=hpsd3.Data;
　　% Instantiate a PSD data object and store the three different
　　% estimates since they all share the same frequency information.
　　hpsd=dspdata.psd(［Pxx1, Pxx2, Pxx3］,W,'Fs',Fs);
　　plot(hpsd);
　　legend('Hamming','kaiser','Chebyshev');
　　end
　　3 结论
　　用产生白噪声加滤波器的方法比较容易仿真产生1/f噪声，并且仿真得到的1/f噪声较真实。
　　
　　参考文献
　　［1］Tae-Hoon Lee, Gyuseong Cho Chap1. Monte Carlo based time-domain Hspice noise simulation for CSA-CRRC cricuit［J］. Nuclear Instruments and Methods in Physics Research Section A, 2003,505:328-333.
　　［2］Terry,S.C. Blalock, B.J. Rochelle, J.M. Ericson, M.N. Caylor, S.D. Time-domain noise analysis of linear time-Invariant and linear time-variant systems using MATLAB and HSPICE［C］. Nuclear Science, IEEE Transactions on 2005,52:805-812.
　　［3］王爱萍,王惠南.基于小波分析的1/f噪声降噪 ［J］.数据采集与处理,2006.
　　［4］文 俊,陈良栋.1/f噪声源以及MOS器件中的1/f噪声［J］.电子器件,1996.