摘 要：本文介紹了1/f噪聲的含義和原理，依照這些原理，通過白噪聲加簡化得到的有理1/f數(shù)字濾波器的方法，使用Matlab編程，仿真產(chǎn)生了真實1/f噪聲。并根據(jù)常見的計算功率譜密度的算法，驗證了仿真產(chǎn)生的1/f噪聲的性能。

關(guān)鍵詞:噪聲；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

近年來，隨著在物理、化學、生物醫(yī)學等基礎(chǔ)學科中對自然現(xiàn)象及規(guī)律的深入探索，在通信、測量、導(dǎo)航、醫(yī)療儀器等工業(yè)領(lǐng)域中人們必須測量及接受越來越小的電信號。測量或者接受裝置中的內(nèi)部噪聲影響越來越為明顯，這種噪聲已經(jīng)成為制約電子儀器和裝備進一步提高信號檢測靈敏度及改進接受信號質(zhì)量的關(guān)鍵因素。仿真產(chǎn)生噪聲信號可以幫助研究人員研究克服噪聲的方法，幫助技術(shù)人員研制出高抗干擾的儀器設(shè)備或指出改進產(chǎn)品的方向。

1/f噪聲廣泛存在于電子元器件和電子線路中，從1/f噪聲測量分析中可以得到很多有用的信息。本文主要討論1/f噪聲的特點和仿真產(chǎn)生的方法。

Matlab的信號處理工具箱是信號算法文件的集合，它處理的對象是信號和系統(tǒng)，信號處理工具箱提供了很多命令行函數(shù)來進行譜分析，包括經(jīng)典的無界技術(shù)以及特征值技術(shù)。另外，也添加了很多對象，用以提高可用性。因此我們使用該工具箱來仿真產(chǎn)生1/f噪聲源，并且用其來分析驗證仿真產(chǎn)生出的噪聲可信度。

1 1/f噪聲的特性

1.1 1/f噪聲（低頻率噪聲）的特性

電子器件中的1/f 噪聲具有兩個基本特性：

1)1/f 噪聲在一個相當寬的頻帶內(nèi)，功率譜密度與頻率f 成反比， 且頻帶上、下限都是有限值，上限頻率視1/f 噪聲與背景白噪聲的相對大小而定；下限頻率接近直流時功率密度仍呈很好的1/f 特性。

2)1/f 噪聲電壓和電流的功率譜密度近似與流經(jīng)器件的電流的平方成正比，這意味著1/f噪聲起源于電阻的漲落。設(shè)流過電阻R的電流保持恒定，電阻存在著漲落，引起電阻兩端電壓漲落，則有:

δV=1*δR(1-1)

Sv(f)=I2SR(f)I2(1-2)

假設(shè)功率譜密度是Pnf=Sf/f，相應(yīng)的電壓密度是e2nf=Af/f。從頻率f1到f2之間的總噪聲是

v2nf=f2f1Affdf=Aflog(f2f1)(1-3)

所以，1/f譜圖取決于上下的截止頻率，而不是絕對帶寬。

1.2 1/f噪聲的自相關(guān)函數(shù)

LTV(線性時不變)系統(tǒng)的輸出可以用卷積積分來表示：

(1-4)

式中t是輸出的觀察時間。

通常我們更關(guān)心的是計算系統(tǒng)輸出的均方根值，假設(shè)輸入是由一個靜態(tài)噪聲源驅(qū)動的，為了得到這個結(jié)果，第一步是計算輸出均方根值

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次抽樣函數(shù)。

式(1-7)中間括號的內(nèi)容是隨機輸入噪聲信號的自相關(guān)函數(shù)，所以可以簡化成為

x-2∞-∞∞-∞h(u)h(v)m(t-u)m(t-v)［Rnn(u-v)］dudv（1-8）

上面這個式子指出了自相關(guān)函數(shù)在時域噪聲分析中的重要性。在噪聲處理中最常見的兩種噪聲是白噪聲和1/f噪聲的自相關(guān)函數(shù)。理想的1/f噪聲是一個多樣的過程，它的自相關(guān)函數(shù)無法用閉合的形式來簡單的表達出來。實際應(yīng)用中進行了一些近似。

白噪聲的雙邊功率譜密度(double-sidedPSD)是

v2Δf=N02（1-9）

N0是單邊噪聲功率譜密度，單位是V2/HZ。自相關(guān)函數(shù)就是簡單的功率譜密度的富利葉反變換。所以，白噪聲的自相關(guān)函數(shù)由沖擊響應(yīng)函數(shù)給出

Rnn(τ)=N02δ(τ)(1-10)

1/f噪聲，可以看成是白噪聲通過一個噪聲濾波器得到的結(jié)果。白噪聲驅(qū)動的通過線性時間濾波器的輸出的自相關(guān)函數(shù)是

Ryy(τ)=N02∞-∞h(λ)h(τ+λ)dλ(1-11)

h(t)是成型濾波器的沖擊響應(yīng)。理想的1/f成型噪聲濾波器得沖擊響應(yīng)是

h(t)2fct，t＞0(1-11)

fc是轉(zhuǎn)角頻率。1/f噪聲的自相關(guān)函數(shù)從下式計算得到：

Ryy(τ)=N0fc∞01λ1λ+τdλ(1-13)

可以得到

Ryy=2N0fc1n［λ+λ+τ］∞0(1-14)

上面這個式子形式是非常簡單的，可惜的是，式(1-14)是發(fā)散的，所以并不可能去精確的計算得到1/f噪聲的自相關(guān)函數(shù)。我們可以簡化得到有理形式的濾波器（或者一串濾波器），解決方法是定義一個總的操作時間，并使用恰當?shù)慕?，?1-14)可以轉(zhuǎn)變?yōu)?/p>

Ryy(top，τ)=2N0fc1n［λ+λ+τ］top-τδ(1-15)

假設(shè)topτδ，可以得到

Ryy(top，τ)≌N0fc［1n(4)+1n(topfc）-1n(τfc) ］ （1-16）

注意，top和δ，是等同于時域內(nèi)的兩個頻率限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噪聲仿真的實際應(yīng)用中，式(1-14)比式(1-19)更實用。

1.3 1/f噪聲的產(chǎn)生方法

產(chǎn)生1/f噪聲的最直觀的方法就是將白噪聲通過噪聲濾波器過濾之后得到。但是，理想的1/f濾波器并不是有理的，所以需要開發(fā)一個可以近似得到1/f噪聲的濾波器。需解決的主要問題是求這一近似的濾波器的算法。

使用下圖所示簡化的噪聲濾波器。

圖1 仿真1/f噪聲成型濾波器的示意圖

一系列的hn(t)進行疊加，可以得到一系列的濾波器，產(chǎn)生了完整的1/f噪聲的輸出，傳輸函數(shù)如下給出

H(s)=1+∑N+1n=010-0.5n22πfcs+2π10-nfc(1-19)

N定義了總頻率帶寬的10倍頻數(shù)量，fc是轉(zhuǎn)角頻率。1/f噪聲的數(shù)字化產(chǎn)生需要將模擬濾波器變成數(shù)字化濾波器，z域傳輸函數(shù)很容易計算得到，計算的結(jié)果如下式：

H(z)=1+∑Nn=010-0.5n2πfcTs1-z-1exp(-2π10-nfcTs)(1-20)

現(xiàn)在可以使用Matlab內(nèi)置的filter函數(shù)來進行處理了。在仿真產(chǎn)生1/f噪聲時，比較合理N的值大約為10。這意味著一個帶寬為1MHz的系統(tǒng)，最低頻率在10-4左右。這樣，總的仿真時間長度需要在s。一般來說104s，仿真N個10倍頻段需要大約10N個點。如果N>6，考慮到對計算機內(nèi)存和資源的要求，是不現(xiàn)實的。在實際應(yīng)用中，通常選定采樣頻率為107，模擬時間長度為10-3?？紤]到所需要的精度大概需要104-105個點。為了解決這個問題我們定義一個參量Neff

Neff=1+log10(tsimTs)(1-21)

tsim是總的仿真時長。然后整個數(shù)字濾波器被分為兩個部分，如圖2所示，所有的nNeff子濾波器，按照式(1-20)來計算，所有的n＞Neff部分，輸出在整個仿真過程中幾乎都是常量。意味著高階的子濾波器可以進行簡單的將數(shù)值（偏移量）相加。每一部分所得均方根等于穩(wěn)態(tài)時子濾波器的輸出均方根。

圖2 數(shù)字濾波器可以被Neff分為兩個部分

1.4 產(chǎn)生白噪聲并計算其功率譜密度

1.4.1 使用matlab仿真產(chǎn)生白噪聲

如上節(jié)分析，因為仿真產(chǎn)生1/f噪聲的方法是白噪聲通過一個濾波器，所以先需要仿真得到一個白噪聲。(程序較簡單略)

在Matlab中產(chǎn)生白噪聲，可以使用randn函數(shù)的得到一個具有高斯正態(tài)分布的白噪聲

noise=rms×randn(1，npts)

產(chǎn)生得到一個均值為0，均方根值為rms，總點數(shù)為npts的噪聲，均方根值由白噪聲底階和采樣頻率共同決定，如下式得到：

rms=N02Ts(1-22)

N0是白噪聲的單邊功率譜密度，Ts是采樣頻率。采樣頻率最好取為10的因子并且略大于系統(tǒng)帶寬，這樣在整個帶寬內(nèi)白噪聲都有一個平坦的功率譜。

程序輸入輸出量說明：

function noise=shot_noise(Ts，N0，npts)

Ts采樣時間，fs=1/Ts采樣頻率，N0白噪聲的單邊功率譜密度，npts仿真所取的點，noise為所得到的噪聲數(shù)值序列

圖3 仿真得到的白噪聲時域波形

1.4.2 計算功率譜密度的三種不同的方法

函數(shù)中使用direct關(guān)鍵詞則使用了直接法進行計算，又稱周期圖法，它是把隨機序列x(n)的N個觀測數(shù)據(jù)視為一能量有限的序列，直接計算x(n)的離散傅立葉變換，得X(k)，然后再取其幅值的平方，并除以N，作為序列x(n)真實功率譜的估計。

使用indirect關(guān)鍵詞，使用間接法先由序列x(n)估計出自相關(guān)函數(shù)R(n)，然后對R(n)進行傅立葉變換，得到x(n)的功率譜估計。

其中complex關(guān)鍵詞，綜合方法會將加入三種窗函數(shù)(Hamming，Kaiser，Chebyshev)的圖譜估計直接顯示在一幅圖上，比較直觀和明顯。

1.4.3 白噪聲的功率譜密度

為了驗證產(chǎn)生的白噪聲的性能，計算它的PSD進行驗證

function psd_s(Ts，N0，npts，method)

% PSD of shot noise

除了method外與產(chǎn)生程序的變量意義相同，method為產(chǎn)生功率譜密度算法參數(shù)。

圖4 仿真得到的白噪聲功率譜

1.5 仿真產(chǎn)生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倍頻數(shù)輛

% tsim:總仿真時長

% fs:采樣周期

% fc:轉(zhuǎn)角頻率

具體源程序代碼見2.1.1節(jié)

圖5 仿真得到的1/f噪聲時域波形

tsim，仿真持續(xù)時間，Ts采樣時間，N10倍頻段，fs采樣頻率，npts總共仿真的點數(shù)。

得到了一個雜亂無章的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節(jié)。

除了method外與產(chǎn)生程序的變量意義相同，method為產(chǎn)生功率譜密度算法參數(shù)。

direct用直接法進行計算，直接法又稱周期圖法，它是把隨機序列x(n)的N個觀測數(shù)據(jù)視為一能量有限的序列，直接計算x(n)的離散傅立葉變換，得X(k)，然后再取其幅值的平方，并除以N，作為序列x(n)真實功率譜的估計。

ndirect，間接法由序列x(n)估計出自相關(guān)函數(shù)R(n)，對R(n)進行傅立葉變換，得到x(n)的功率譜估計。

其中complex方法，會將加入三種窗函數(shù)的圖譜估計直接顯示在一幅圖上。

圖6 仿真得到的1/f噪聲功率譜密度

可以看到，PSD近似與lnf成一個反比的關(guān)系，隨著f的增大而減小。對圖形進行處理，橫軸以ln（1/f）為坐標，下降的關(guān)系是成直線的，在低頻率范圍內(nèi)出現(xiàn)的扭曲，是濾波器傳遞函數(shù)在N＞Neff時的近似造成的。

圖7 對數(shù)坐標-1/f功率譜密度

取對數(shù)坐標作圖(見圖7)，舍棄超過轉(zhuǎn)角頻率的點，進行線性擬合得到：b=-1.4486 a=-2.2333r=-3.9922

2 仿真過程的實現(xiàn)

2.1 仿真產(chǎn)生1/f噪聲

2.1.1 仿真產(chǎn)生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，F(xiàn)s);

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，F(xiàn)s，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'，F(xiàn)s，'SpectrumType'，'twosided'，'CenterDC'，true);

psd(h，xn，hopts)

case{'kaiser'}

h=spectrum.welch('hamming'，64);

hpsd1=psd(h，xn，'Fs'，F(xiàn)s);

W=hpsd1.Frequencies;

h.WindowName='Kaiser';

hpsd2=psd(h，xn，'Fs'，F(xiàn)s);

Pxx2=hpsd2.Data;

hpsd=dspdata.psd(Pxx2，W，'Fs'，F(xiàn)s);

plot(hpsd);

legend('kaiser');

case{'chebyshev'}

h=spectrum.welch('hamming'，64);

hpsd1=psd(h，xn，'Fs'，F(xiàn)s);

W=hpsd1.Frequencies;

h.WindowName='Chebyshev';

hpsd3=psd(h，xn，'Fs'，F(xiàn)s);

Pxx3=hpsd3.Data;

hpsd=dspdata.psd(Pxx3，W，'Fs'，F(xiàn)s);

plot(hpsd);

legend('Chebyshev');

case{'complex'}

h= spectrum.welch('hamming'，64);

% Create three power spectral density estimates.

hpsd1=psd(h，xn，'Fs'，F(xiàn)s);

Pxx1=hpsd1.Data;

W=hpsd1.Frequencies;

h.WindowName='Kaiser';

hpsd2=psd(h，xn，'Fs'，F(xiàn)s);

Pxx2=hpsd2.Data;

h.WindowName='Chebyshev';

hpsd3=psd(h，xn，'Fs'，F(xiàn)s);

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'，F(xiàn)s);

plot(hpsd);

legend('Hamming'，'kaiser'，'Chebyshev');

end

3 結(jié)論

用產(chǎn)生白噪聲加濾波器的方法比較容易仿真產(chǎn)生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］.數(shù)據(jù)采集與處理，2006.

［4］文 俊，陳良棟.1/f噪聲源以及MOS器件中的1/f噪聲［J］.電子器件，1996.