１／ｆ噪聲的仿真產生和驗證

摘 要：本文介紹了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.