Generation of 1/f noise using Matlab.

Dear Colleagues, I have been trying to generate the 1/f noise, where f means frequency. I would appreciate any help and guidance. Kind regards,

Massilon

Answer 1

[Star Strider]

Probably the easiest way is to create a FIR filter that has a ‘1/f’ passband, then filter random noise through it:

fv = linspace(0, 1, 20);                                % Normalised Frequencies
a = 1./(1 + fv*2);                                      % Amplitudes Of ‘1/f’
b = firls(42, fv, a);                                   % Filter Numerator Coefficients

figure(1)
freqz(b, 1, 2^17)                                       % Filter Bode Plot

N = 1E+6;
ns = rand(1, N);
invfn = filtfilt(b, 1, ns);                             % Create ‘1/f’ Noise

figure(2)
plot([0:N-1], invfn)                                    % Plot Noise In Time Domain
grid

FTn = fft(invfn-mean(invfn))/N;                         % Fourier Transform

figure(3)
plot([0:N/2], abs(FTn(1:N/2+1))*2)                      % Plot Fourier Transform Of Noise
grid

It uses the firls function to design a FIR filter that closely matches the ‘1/f’ passband. See the documentation on the various functions to get the result you want.

Note: The filter is normalised on the open interval (0,1), corresponding to (0,Fn) where ‘Fn’ is the Nyquist frequency, or half your sampling frequency. It should work for any sampling frequency that you want to use with it.

This should get you started. Experiment to get the result you want.

[Massilon Toniolo da Silva]

I have got almost a perfect 1/x plot. Excellent, superb! Great job! Many thanks, Massilon.

[Star Strider]

As always, my pleasure!

Simulating ‘1/f’ noise is an interesting problem that I had never before considered.

[Hour Elen Sayed]

Could I ask what does those normalized frequencies represent ?

[Star Strider]

The frequencies are normalised with respect to the Nyquist frequency. The normalised frequencies represent radian frequencies in the range of radians/time unit.

[Antonio D'Amico]

What if I want to add a white noise component, so to play with the corner frequency?

[Star Strider]

Use the randn function instead of rand:

ns = randn(1, N);

Make appropriate changes to get the result you want.

[Antonio D'Amico]

Hello, thank you for your answer. If I understand the script correctly, it applies a 1/f (approximation) roll-off factor to the noise, whether it is uniformilly distributed (rand) or gaussian (randn). However what I would like to achieve is something like (From Wikimedia Commons)

I hope I was clearer

[Antonio D'Amico]

Ok, I think I got it, something like this could work

fv = linspace(0, 1, 20);                                % Normalised Frequencies
a = zeros(1,20);
a(1:10) = 1./(1 + fv(1:10)*2);                          % Amplitudes Of  1/fv until 0.5
a(11:20) = a(10);                                       % after 0.5 it gets flat
b = firls(42, fv, a);                                   % Filter Numerator Coefficients
figure(1)
freqz(b, 1, 2^17)                                       % Filter Bode Plot
N = 1E+6;
ns = randn(1, N);
invfn = filtfilt(b, 1, ns);                             % Create ‘1/f’ Noise
figure(2)
plot([0:N-1], invfn)                                    % Plot Noise In Time Domain
grid
FTn = fft(invfn-mean(invfn))/N;                         % Fourier Transform
figure(3)
plot([0:N/2], abs(FTn(1:N/2+1))*2)                      % Plot Fourier Transform Of Noise
grid

Answer 2

[Ali Mostafa]

f=0:1/fs:1-1/fs;S=1./sqrt(f); S(end/2+2:end)=fliplr(S(2:end/2)); S=S.*exp(j*2*pi*rand(size(f))); plot(abs(S)) S(1)=0;figure;plot(real(ifft(S)))

[Massimo Ciacci]

Quite ingenious to put the randomness in the phase, and this way the amplitude profile is exact, without the need to average out a lot of noise realizations. Thumbs up!

[XIAOHUA HUA]

Great, thank you very much for sharing this.

Answer 3

[James]

Hi were does 1./(1 + fv*2) come from?

[Star Strider]

It keeps the amplitude vector from becoming infinite at the origin. With that amplitude calculation, it is 1 at the origin. The factor of 2 (that can be anything that works in a particular application), allows the amplitude to decrease differently than simply . If the factor is less than 1, the decrease is slower, if greater than 1, faster.

[James]

is there any paper or book I could look at to undestand that a bit more, or is this based on your own experience/skill?

Thank you very much for your response!

[Star Strider]

It’s entirely my own experience. I remember learning about noise in graduate school, in the context of biomedical instrumentation. I’m certain there must be more recent discussions of it, however I have no specific references.