跪求生成粉红噪声的MATLAB程序

　MATLAB里面没有专门产生这种噪声的命令，不过你可以先用rand, randn这陪颤则样的产生随机白噪声，然后用filter来芦棚滤一下。就能得到粉噪声了。

% function [pn, theta] = phase_noise(num_samp, f0, dbc_per_hz, num_taps)

%

% This function creates noise with a 1/f spectrum. The noise is then

% phase modulated so that it can be mixed with a signal to simulate

% phase noise in the original signal. The noise is specified in power

% per hertz at frequency f0, relative to the power in the carrier.

%

% References:

% N. J. Kasdin, "Discrete Simulation of Colored Noise and Stochastic

% Processes and 1/f^a Power Law Noise Generation," _Proceedings of

% the IEEE_, May, 1995.

% Roger L. Freeman, _Reference Manual for Telecommunications

% Engineering_.

% M. Schroeder, _Fractals, Chaos, and Power Laws_.

%

% Input/Output parameters:

% num_samp desired number of output samples

% f0 reference frequency (must be in Hz.)

% dbc_per_hz power per hertz relative to carrier at ref. freq.

% num_taps number of filter taps in AR 1/f filter

% (optionaldefault = 100)

%

% pn phase-modulated 1/f process

% theta 1/f process (before phase modulation)

%

% Jeff Schenck 11/21/95

%

% 1/f noise is produced by passing white noise through a filter. The

% resulting spectrum has the form

%

% Sx(w) = g^2 / w, (pretend that w is an omega)

%

% where g is the gain applied to the white noise before filtering. If P

% is the desired power in the 1 Hz. band at w0 = 2pi*f0/fs, and W is the

% 1 Hz. bandwidth in radians, we can write

%

% P = (1/2pi) Sx(w0) W

% = (1/2pi) g^2/w0 (2pi*1/fs)

% = g^2 /洞慧 2pi*f0

%

% =>g = sqrt(2pi*f0*P).

%

% Notice that the result is *independent* of fs!! Look at it this way:

% if the sampling rate is doubled for a given spectrum, the new w0 is

% half the old w0. For 1/f noise, this means that Sx(w0_new) =

% 2*Sx(w0_old). But a 1 Hz. band is half as large (in radians) than it

% was previously, so the product P is the same, and fs drops out of the

% picture.

%

% The independence with respect to fs is also an indication of the

% fractal nature of pink noise.

%

% Note that the phase-modulated noise is itself 1/f if the narrowband

% assumption is valid.

function [pn, theta] = phase_noise(num_samp, f0, dbc_per_hz, num_taps)

% Check input.

if dbc_per_hz >= 0

error('Power per Hz. must be negative.')

elseif f0 <= 0

error('Reference frequency must be positive.')

end

if nargin <4

num_taps = 100

end

% Generate white noise. Apply gain for desired dBc/Hz. Warn user

% if gain is too large (gain thresholds have been chosen somewhat

% arbitrarily -- needs work).

gain = sqrt(2*pi * f0 * 10^(dbc_per_hz/10))

wn = gain * randn(1,num_samp)

fprintf('Gain applied to white noise = %f.\n', gain)

if gain >= 1

fprintf('WARNING: Narrowband approximation no longer valid.\n')

elseif gain >= .5

fprintf('WARNING: Narrowband approximation on the verge of collapse.\n')

end

% Generate 1/f AR filter and apply to white noise to produce 1/f

% noise.

a = zeros(1,num_taps)

a(1) = 1

for ii = 2:num_taps

a(ii) = (ii - 2.5) * a(ii-1) / (ii-1)

end

theta = filter(1,a,wn)

% Phase modulate.

pn = exp(i*theta)

return

宇宙中的“1/f 波动(1/f fluctuation ，1/f 噪声)”具有维系大环境平衡、和谐的重要作用，在国外，“1/f 波动”引洞禅起了科学家们的极大重视。自然界中不存在绝对静止的事物，一切事物必定随时间或空间的推移而发生变化，我们将某个物理量在宏观平均值附近的随机变化称为“波动”。波动现象随处可见：有直观的，如风速的快慢变化、声音的节奏和强弱变化；也有抽象的， 如人的情绪变化。自然界存在着许多波动，可以按功率谱密度与频率的对应关系对其进行分类。有三种典型噪声的波动特性。白噪声是一种完全无规律的令人烦躁不安的噪声，该噪声的功率谱密度平行于横轴，是与频率无关的量，我们称之为 1/f0 波动；布朗噪声是一种相关性很强，使人感到单调乏味的噪声，该噪声的功率谱纳岁尘密度与 f2 成反比，我们称之为 1/f2 波动；介于上述两种形式之间的噪声是一种在局部呈无序状态，而在宏观上具有一定相关性的噪声，是一种使人感到舒服的波动。由于该噪声的功率谱密度与频率在是成反比的，我们称之为 1/f 波动。事实证明：“1/f 波动”与人在安静时的α脑波及心拍周期等生物体信号的变化节奏相吻合并与人的情感、感觉有着密切联系，使人能感到舒适。人的感觉千差万别造成对美的感受各有不同，我们缺少一个雀键客观评估的标准。而“1/f 波动”理论的提出使我们有了这样一个准则。

很简单。一般是先产生白噪声信号，傅里叶变换到功率频谱，对白噪声功率谱乘以pink噪声的功率谱函数，就是频卜镇域处理，然后再反向傅里叶变换到时型锋粗域基此波形。粉色噪声的频谱是1/f。在频域比较好处理。记住频域信号的相乘等效于时域波形的互相关，这个是信号处理的基本套路。

---