Generate pinknoise with same power as audio

조회 수: 2 (최근 30일)
L . 2023년 8월 25일
답변: Mathieu NOE . 2023년 8월 28일
I need to read an audio file and generate a pinknoise with the same power.
I tought using generatepinknoise(), but the built-in function does not accept a power specification.
Any ideas?
I was trying something based on, but it dosen't sound anything like pinknoise.
% Generate pinknoise with same power as sound
% get sound (2 channels)
[y,Fs] = audioread('S5.wav');
A1 = abs(Y(:,1)); P1 = mean(A1.^2);
A2 = abs(Y(:,2)); P2 = mean(A2.^2);
N = length(Y); invfnorm = 1./[1:N]';
Anew1 = sqrt(P(:,1)*invfnorm/sum(invfnorm)); Anew2 = sqrt(P(:,2)*invfnorm/sum(invfnorm));
noise = ifft([Anew1 .* exp(1i*angle(Y(:,1))) Anew1 .* exp(1i*angle(Y(:,1))) ] ,'symmetric');

채택된 답변

Mathieu NOE
Mathieu NOE 2023년 8월 28일
this would be my suggestion . You don't need to do a fft to get the power of a signal (you can compute it directly in time domain too)
here I tested it on my own audio signal (piano.wav)
% Generate pinknoise with same power as sound
% get sound (2 channels)
[y,Fs] = audioread('S5.wav'); % your file
% [y,Fs] = audioread('piano.wav'); % mine
[samples,channels] = size(y);
y_power = mean(y.^2,1); % power of signals
pn = mypinknoise(samples);
pn_power = mean(pn.^2); % yes it's 1 because normalized std as explained in the subfunction
% let's assume we want same pink noise power as audio channel of largest
% power
pow_factor = max(y_power)/pn_power;
% now the pink noise signal must be multiplied by the square root of the
% power factor
pn = pn*sqrt(pow_factor);
pn_power = mean(pn.^2); % double check ; now we have same noise power as the audio signal
% concatenate signals for replay
out = [y; pn(:)*ones(1,channels)];
function y = mypinknoise(N)
% function: y = pinknoise(N)
% N - number of samples to be returned in a row vector
% y - a row vector of pink (flicker) noise samples
% The function generates a sequence of pink (flicker) noise samples.
% In terms of power at a constant bandwidth, pink noise falls off at 3 dB/oct, i.e. 10 dB/dec.
% difine the length of the vector and
% ensure that the M is even, this will simplify the processing
if rem(N, 2)
M = N+1;
M = N;
% generate white noise
x = randn(1, M);
X = fft(x);
% prepare a vector with frequency indexes
NumUniquePts = M/2 + 1; % number of the unique fft points
n = 1:NumUniquePts; % vector with frequency indexes
% manipulate the left half of the spectrum so the PSD
% is proportional to the frequency by a factor of 1/f,
% i.e. the amplitudes are proportional to 1/sqrt(f)
X = X(1:NumUniquePts);
X = X./sqrt(n);
% prepare the right half of the spectrum - a conjugate copy of the left one,
% except the DC component and the Nyquist component - they are unique
% and reconstruct the whole spectrum
X = [X conj(X(end-1:-1:2))];
y = real(ifft(X));
% ensure that the length of y is N
y = y(1, 1:N);
% ensure unity standard deviation and zero mean value
y = y - mean(y);
y = y/std(y, 1);

추가 답변 (0개)


Help CenterFile Exchange에서 Measurements and Spatial Audio에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by