Generation of 1/f noise using Matlab.

136 ビュー (過去 30 日間)
Massilon Toniolo da Silva
Massilon Toniolo da Silva 2017 年 6 月 14 日
コメント済み: Star Strider 2023 年 10 月 3 日
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

採用された回答

Star Strider
Star Strider 2017 年 6 月 14 日
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.
  8 件のコメント
Antonio D'Amico
Antonio D'Amico 2020 年 8 月 26 日
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)
(From Wikimedia Commons)
I hope I was clearer
Antonio D'Amico
Antonio D'Amico 2020 年 8 月 26 日
編集済み: Antonio D'Amico 2020 年 8 月 26 日
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

サインインしてコメントする。

その他の回答 (2 件)

Ali Mostafa
Ali Mostafa 2018 年 6 月 11 日
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)))
  2 件のコメント
Massimo Ciacci
Massimo Ciacci 2019 年 8 月 10 日
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
XIAOHUA HUA 2020 年 3 月 11 日
Great, thank you very much for sharing this.

サインインしてコメントする。


James
James 2023 年 10 月 3 日
Hi were does 1./(1 + fv*2) come from?
  3 件のコメント
James
James 2023 年 10 月 3 日
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
Star Strider 2023 年 10 月 3 日
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.

サインインしてコメントする。

カテゴリ

Help Center および File ExchangeDigital Filter Analysis についてさらに検索

タグ

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by