From 2e32b8d1f84e484c2635a2f836e92bc08c6cbddd Mon Sep 17 00:00:00 2001 From: wizardforcel <562826179@qq.com> Date: Wed, 11 Apr 2018 11:08:49 +0800 Subject: [PATCH] 8. --- 8.md | 88 +++++++++++++++++++++++++++++++++++++++++++++++++++- img/8-6.png | Bin 0 -> 6568 bytes 2 files changed, 87 insertions(+), 1 deletion(-) create mode 100644 img/8-6.png diff --git a/8.md b/8.md index 28ffa42..81545fe 100644 --- a/8.md +++ b/8.md @@ -20,7 +20,7 @@ + 一些物理量的重尾分布：例如，在冷冻的水中，晶体尺寸的分布是幂律的。 -+ 呈现粉红噪音的时间变化。复合信号可以分解为它们的频率分量。在粉红噪声中，低频分量比高频分量功率更大。具体而言，频率`f`处的功率与`1 / f`成正比。 ++ 呈现粉红噪声的时间变化。复合信号可以分解为它们的频率分量。在粉红噪声中，低频分量比高频分量功率更大。具体而言，频率`f`处的功率与`1 / f`成正比。 临界系统通常不稳定。例如，为了使水保持部分冷冻状态，需要主动控制温度。如果系统接近临界温度，则小型偏差倾向于将系统从一个相位移到另一个相位。 @@ -369,3 +369,89 @@ slope = params[0] 严格来说，值为 3 的分形维度与 2 不可区分，但考虑到其他值的结果，线的表观曲率和图案的外观，似乎它也是分形的。 本章结尾的练习之一，要求您使用不同的`n`和`level`值，再次运行此分析，来查看估计的维度是否一致。 + +## 8.6 频谱密度 + +提出沙堆模型的原始论文的标题是《自组织临界：1/f 噪声的解释》（Self-Organized Criticality: An Explanation of 1/f Noise）。 正如小标题所述的那样，Bak，Tang 和 Wiesenfeld 正试图解释为什么许多自然和工程系统表现出 1/f 噪声，这也被称为“闪烁噪声”和“粉红噪声”。 + +为了了解粉红噪声，我们必须绕路来了解信号，频谱分析和噪声。 + +信号： + +信号是随时间变化的任何数量。 一个例子是声音，即空气密度的变化。 在本节中，我们将探讨雪崩持续时间和大小在不同时间段内的变化。 + +频谱分析： + +任何信号都可以分解为一组具有不同音量或功率的频率分量。 例如，演奏中央 C 上方的 A 的小提琴的声音，包含频率为 440 Hz 的主要分量，但它也包含较低功率分量，例如 880 Hz，1320 Hz 和其他整数倍的基频。 频谱分析是寻找构成信号的成分和它们的功率的过程，称为其频谱。 + +噪声： + +在通常的用法中，“噪声”通常是一种不需要的声音，但在信号处理的情况下，它是一个包含许多频率成分的信号。 + +噪声有很多种。 例如，“白噪声”是一个信号，它在很宽的频率范围内拥有相同功率的成分。 + +其他种类的噪声在频率和功率之间有不同的关系。 在“红噪声”中，频率为`f`的功率为`1 / f ** 2`，我们可以这样写： + +``` +P(f) = 1 / f ** 2 +``` + +我们可以把指数 2 换成`β`来扩展它： + +``` +P(f) = 1 / f ** β +``` + +当`β= 0`时，该等式描述白噪声; 当`β= 2`时，它描述红噪声。 当参数接近 1 时，我们将结果称为`1 / f`噪声。 更一般来说，任何介于 0 和 2 之间的噪声称为“粉红”，因为它介于白色和红色之间。 + +那么这如何适用于沙堆模型呢？ 假设每次细胞倒塌时，它会发出声音。 如果我们在运行中记录了沙堆模型，它会是什么样子？ + +在我的`SandPile`实现运行时，它会记录在每个时间步骤中，倒塌的细胞数量，并将结果记录在名为`toppled_seq`的列表中。 例如，在第？节中运行模型之后，我们可以提取产生的信号： + +```py +signal = pile2.toppled_seq +``` + +为了计算信号的频谱（同样，这是它包含的频率和它们的功率），我们可以使用快速傅立叶变换（FFT）。 + +唯一的问题是噪声信号的频谱往往是嘈杂的。 但是，我们可以通过将一个长信号分成多个段，计算每个段的 FFT，然后计算每个频率的平均功率来使其平滑。 + +该算法的一个版本被称为“韦尔奇方法”，SciPy 提供了一个实现。 我们可以像这样使用它： + +```py + +from scipy.signal import welch + +nperseg = 2048 +freqs, spectrum = welch(signal, nperseg=nperseg, fs=nperseg) +``` + +`nperseg`是信号分解成的片段长度。 对于较长的片段，我们可以获得更多的频率成分，但由于平均的片段数较少，结果更加嘈杂。 + +`fs`是“采样频率”，即每单位时间的信号中的数据点数。 通过设置`fs = nperseg`，我们可以得到从 0 到`nperseg / 2`的频率范围，但模型中的时间单位是任意的，所以频率并不意味着什么。 + +返回值，`freqs`和`powers`是 NumPy 数组，包含成分的频率及其相应的功率。 + +如果信号是粉红噪声，我们预计： + +``` +P(f) = 1 / f ** β +``` + +对两边取对数会得到： + +``` + +logP(f) = −β logf +``` + +所以如果我们在双对数刻度上绘制功率与频率的关系，我们预计有一条斜率为`β`的直线。 + +![](img/8-6.png) + +图 8.6：随着时间推移的倒塌细胞的功率频谱，双对数刻度 + +图？显示结果。 对于 10 到 1000 之间的频率（以任意单位），频谱落在一条直线上。 灰线斜率为 -1.58，这对应于由 Bak，Tang 和 Wiesenfeld 报告的参数`β= 1.58`。 + +这个结果证实了沙堆模型产生粉红噪声。 + diff --git a/img/8-6.png b/img/8-6.png new file mode 100644 index 0000000000000000000000000000000000000000..e9c2870203b5c388d71eb4623df7a33052718202 GIT binary patch literal 6568 zcmZ`;XHXMNw+^Bd6_FxRLI){=LTCxSgeJXL0Z9UePUuBYihw}qNC~}4l_phqX$ieY zLWLXlQ7Bd^{$*rK_u}cVrg%u4sPo zSA2GTBBreVOS7(1s`9I7D^z?=No{0gB-AF#0Qn&$w?;%nq^zurgM&j0nWXUYy#_2I zF{f_t_u&h>I9c;hW$X95ySsLFcG6}ca+cxIm~%}o}uX% z@MstP*;z(Lrl6nzS@-tj=3qpg z`q~Bn`PJ`s{{mYLRVdN*Ey1}PK<=hy<^cduzWTQimCBVM0RTFGb*O@&?~k3AQ8~+f zM%sv{{S}O<0k^nQFzGH9b(|Kh+x}m%+L6q&(G|!(`WYHZR10(KjKD)S^3Rs+@^9VV zU$o!4<}>1+j}Yij6zJc(*zw0(4UXot`&U~LSkKgabK}7qS=8N*6X9G}xjoiCOYFhT z@!mfQFIV*Xtjqc885gu`_r+zF$O)|bPS%vwA5`MT>S)`!q)*zDVq|ui!=As<`oS9^ z*y!~F`JMXO>EfFV;x$wEZ~#HdS00P7-6=6Cej{JZz2 zI1H^+(nQml#qcR%)k`mzv45lbOsC6rS zBq`yoxpV%L$32my67$dVjU7HkqMh+VtqLRrm`bL0~$Q=@gUHGYZmQucr}&%keFLHDLKByn>Y&<}0Kg&_ijyxk!#EH4SY zNO|vzS~H@V|4aBt9NVg^9P>ytJ~abAhdrAE(HIczL;pHXKM4``u^fz z=Zel`vf=0j_ zZUNmCRK`OEIgUslWqejP++}u7xb2)j*-1smt*HB3w)NImlP+`5u!esK4=I+>sWg-B zbdKZqnzfVa`$%U=lKnMpKFo3G=$?WVI^${T`hKSvCl%dAh0rAkv=}$s7R)c}Ect85 z7$mTyCII5K#DXzIhsyV58F1sqUseHU;GpDUU!}4s`lEb z@mo6C+2)_sARt;l)?lQ^!xNrohdtxU1?na?hT#~Kv!y1YSXUVpJhoua$*qSlmL-V4 zTQnN9diZHIG9gs+YaDt@w`?r)OR6mz_N%I`35-H$xJzEz+R{$jnP6Arr1V;l5-3`a zKku}I&-nL%XR--LE3UsVoDrHO-F$+@vCB#_2dt9F(caps$g1xt<;RJO%W@_tTUcLN{f z6{K@8+m}}LZ4^wgd>aS#Ufvt>)oKwKtZw^x)CUC?X)XN0ZHXhms(ANEIZn@2P&?)^&t4?~Jv2RMH${PU z>1)8nxKCT8zZFTam26f2F+Y!y1e@eF>5^`H_LH;kvsGdOTPCe@YS)Z(-qKiwY z%n4ZLp^2mV$MP`}FWOYt+?BAcDc3cSJH{bwJXJ^TioRPi^Y@dcn`q^jkzbAa#_(uyI;XYi9M6*6;^d>$eTw#YTV`&d4+`JO``V&e*_f*BVWFc*e@zeWOR|%Q2MdgVK&T7@Nc?aj4aQwLP6R@E6D`{N>E)``u0|L?ovikEfr&E#(zU6f~olP*41E zdi!)4Q{3XF0nJn^qBNx_2k)a|Fb;?Muv84einK zT+fZvSTeBp=huv_nQ2qo8A1N&(t*URa0>^VM}3Xx)7&T5^(nWQ2}drH9vX9cW=8OW zp=X<=jPGHUB0gP+oJZ^>JZ6Oh=L%#3J^UO$=)^PPA|M&jM&z)7WH+XnpHEjp=(Apu z%k4%}H^t71JN#|H}6E3a^#}&(;?Bes~ZPTr{Ved1n-LE zh{GU7qosFLVQVnHBc9*r+HM(>4v_SA61kG(n6ovmzV?3kV#cR^v%x_a5ic=YyWL8t z!E1vQ9eti1;HA3*M24*GmE$R0Yr@q!(5)owk0QwfWk7Lzb9 z=mw6!Car2~AIhoUBc#iV7E7c;_c3QFv#+x4i@dJu)H5ci({!$};dPC{DDDlyWc`T3 zhj^}*z~l*t8dYBUUJLBo&!k}VIvJWv14zD31{47--v^&>Z-GtFTn4lX|Bb{gkj=BT z*Grcu%*mRX9okD1r8TnwwyT!oVZ&+uS(7IaK_JOSyyd#a3aVpNwhwMR(?b*XDn&CH z2AnqfqGHu?;MBV_53-cJ1-39yQTZFZdFhc3i9i&MZl;h@1~R*Puoo>zXOR+(Oq!V) zq*899Ay1V%YZZg)Xuo}W!tM%;C3dGpwODMDS`+m%o7Ga0(|}rHY!)I0f;S4I-QS>Ry|;g$)ou)!k#yqzvq?1gF8z-e+PNN zTliQ*Oomp~$s!M&Pq^Ud$s`VIC@uI5Jhb!ys4g^(j$bd$SDj7jif#R}Xocepz z5}%YNv7+-QxX?YsGV{FMw!t9(7Uz}dW%`XUJcjU~!xrdRQT%!P(v7w>YDMlm|5k&u z=m^Xb|7Mz?>jS`v{2~Lu?~uhuyd{n!g5Q1{1Bj;rfI|K&Zpo|si*5v^hg=drU&(7`WX6};$#%i zP5f!BlJ*&iFN`07$uqAg&%{NGJyQNe|6VWI1am;Z?X?qwijA<%85G^6yaTEE90%X- zRb{)j1qErR3o@J-=Jmcn_Nr5yU1C6#pl7)Z0|&WfqnE$z=K&0yNmM>Zg>Xzj$zt}5 zmc{hm)sYk%j24vm45eeQrH%DVd^#04chNkWtblz8qPE?q%g5vFNWCQcrgS(8en#cG zl>pl>{V3q~w_~Zj5#51{pR-K~P0u?jZxplq%JITYF|e?fWvBsgk8BWgu&O~3KklRJ zMgNcyQe0G`tip~6pF1WwDdOl$+a!r#Eits%SsXS@QNYH7WPAK2!Gm7A$&dqerjDFd zD)gytVdWpvj#Fc8xaKvPs?&nO{q}>WI?vfl?Ziq)Y$DOw)j>O39S;=5>gITPlVLU;|36Gq|4hFvlJ@!<|zlt*^CKi+ljkuT3mu0dlic>Cc*P+wEL zu+~>_d%fY$(8oTuj0)_8j-!|W_dxE?DjPGy@kP~lePQnrP>y2bG=reCS;Y1q^;85# z{cJm~{TNS&A5m4&?73g|1G-sj;`IE1U5eiZ@^f~*$1XIg_P>%^@rYnm3d0(TCIH zaul%Zoje_&l0eAI>4ZIv+;=eyj;3X{$nEzU4&nY(@F8l(in78=*9s}_6SeMK2L|{M z9e*n(v625#%0*EzJsV71_gALcYz02qD*?w*6Pl5Fm%)n$F8#I^TOxcya*`>#@$KWr zZTa!&iR z@IY5i!iFbKOtxsi`=s<1gYr}8Q-f)~9ObQeDqonrw0PF;fvap(=yqB(G3YKgVZ5^S zap)g?tiFLDUV~fMX$C$_{U;}%9%D?Ax8J$pOV#(<^->UyX}S)HX*Dp(E_0Du@lS$U zSu&84YUan+Xb9N*sfD*T>IFNg$A4m}csiYhqNo+0XG@*nf_q}K7)8N~DyAvB9s|dM z7>%-4&F3Sl!Olx+Y%pG0Q$f=zw-Tpur7_!-1kWq6>qA+mKutvdiM+>0tiR3PYXny} zcJ;Q;u&nQLf{ctZqCF}sLW*B;D4#A=0&P}iE~0>4X8)1QqHXq=fpU5xGS8)gY6-?6 z+~@nSj$&`GUE~Ng;4clb9ds{b?4zeBvn;Sh)>yrfNRHjFEL8zx*X@6CU%&#JXXSb3skLbBuR^|HWO`bw4c0c7) zUN<+6@+n|77qe9Nc06Ci3wxY4248PKmd9FxSEeN#ZhFc$%3G~u%oyB#!~^exRPHo4 z*NTciVQBm!!Kls|FYFM35sZDqQXo;~>N>^qgopvAp9onAw(RAFu{27&U=TG?Ux~dw z-@8^v1rz|COA{Y3KlWPAEV$Nk3Amx#@_I0g{)s8OoQktlS!Wz=LyujKD(|E<#XK}z zPNWKP_@YPZI;6frHkZaJ|AF`>0d@omPs{sSnw=90T>v!-NfY>(|3^4mft$ki*nRU#)Qo_u#~s80tkFOAMAUl4kF@kIQV3PP_(4 z5xLURn3zTOVw|YiwGB^a(~^>TBtiKWcDjmRkg+rokYZax42 z2AMHo1uXDwL+fWHr=MKjI7Lnk`{1b6;?N`a>D`kA&W)ZWE$1n`I%}RV3Bx2Gr)%Lu zlB~u#9z_*`LeJ9PuxJ*=v>G$=_FS4d!}x0fb=bEf5y7O~5!1~&iGEWMj|jsVYuobR z=8e73boq-esv18l7LqT%%lF{M6rh5fY!|1&mxRKz2NTIRsSE(Or?ROh{FNs{5{5$D z2mnY(7%V6xDkuao6c&{f6O{zrv_&O_gk(@dj|>