Unexpected power spectrum of a signal with uniform random amplitude values

What should the power spectrum of a signal with uniform random amplitude values be? I expect it to be uniform within the constraints of a finite signal, with no frequency dominating.

However, using the following code I do not obtain the expected quasi-uniform power spectrum. In the attached diagram file, the magnitude ressembles a gamma function. Why is this so?

% Spectrum of signal with uniform random amplitude values
clear all

% signal
n = 2^9; % 2^9 = 512, 2^13 = 8192 % signal length
t = 1:n; % time
y = unifrnd(-1,1,[1,n]); % amplitude, uniform random
% mu = 0;
% sigma = 1;
% y = random('Normal',mu,sigma,[1,n]); % amplitude, normal random

% spectral analysis
nyq = ceil(n/2); % Nyquist frequency
f = fft(y); % Fourier transform
f = f(1:nyq); % half the symmetrical spectrum
m = abs(f); % magnitude
m = m.^2; % power
%m = log(m.^2); % log power
%m = sqrt(m);
p = angle(f); % phase

% sorted values
ys = sort(y,'descend');
ms = sort(m,'descend');
ps = sort(p,'descend');

% histograms
hn = 100;
hx = 1:hn;
hy = histcounts(y,hn);
hm = histcounts(m,hn);
hp = histcounts(p,hn);

% visualization
figure('Name',['Analysis of noise spectra (n = ',num2str(n),')'])

% signal
subplot(3,2,1)
hold on
plot(t,y,'b')
plot(t,ys,'r','LineWidth',2)
hold off
xlim('padded')
ylim('padded')
axis square
box on
title('Signal')
xlabel({'Time','Red line: sorted amplitude values.'})
ylabel('Amplitude')

subplot(3,2,2)
stairs(hy,hx,'b')
xlim('padded')
ylim('padded')
axis square
title('Histogram')
xlabel('Counts')
ylabel('Amplitude Values')


% magnitude/power
subplot(3,2,3)
hold on
plot(t(1:nyq),m,'b')
plot(t(1:nyq),ms,'r','LineWidth',2)
hold off
%set(gca,'XScale','log') 
%set(gca,'YScale','log') 
xlim('padded')
ylim('padded')
axis square
box on
%title('Magnitude')
title('Power (Magnitude^2)')
%title('Log Power (Magnitude^2)')
xlabel('Frequency')
%xlabel('Log Frequency')
ylabel('Amplitude')
%ylabel('Log Amplitude')

subplot(3,2,4)
stairs(hm,hx,'b')
xlim('padded')
ylim('padded')
axis square
title('Histogram')
xlabel('Counts')
ylabel('Amplitude Values')


% phase
subplot(3,2,5)
hold on
plot(t(1:nyq),p,'b')
plot(t(1:nyq),ps,'r','LineWidth',2)
hold off
xlim('padded')
ylim('padded')
axis square
box on
title('Phase')
xlabel('Frequency')
ylabel('Radians')

subplot(3,2,6)
stairs(hp,hx,'b')
xlim('padded')
ylim('padded')
axis square
title('Histogram')
xlabel('Counts')
ylabel('Amplitude Values')

Answer 1

[William Rose]

@Vlad Atanasiu,

You are correct that the power spectrum of a random signal is expected to be approximately flat. However, the magntidues of the spectrum, computed as you have done, with no spectrum averaging, have a standard deviation as large as the mean value. This is explained in textbooks on random signals.

If you do spectrum averaging, you reduce the standard deviation, although you lose spectral resolution. You can demonstrate this using the built-in function pwelch().

You have added a plot of the sorted spectral data to the plot of the magnitude spectrum, at subplot(3,2,3). This red line plot is incorrect, because the y values are plotted versus specific frequencies. But when you sort the magntiudes, the frequencies get all mixed up. So the red plot should not be there.

Power spectral estimates from Gaussian noise have a chi-squared distribution. If you use Gaussian noise instead of uniform, and make a histogram of the spectrum, you will confirm that this is true. I have done it in the past.

Good luck with your work.

[William Rose]

@Vlad Atanasiu,

The code below demonstrates that the spectrum of random noise is similar at all frequencies, and this becomes very clear, if you do spectrum averaging with the Welch periodogram. See the help for pwelch() for an explanation.

N=2^16; x=rand(1,N)-0.5;    %x=random ~U[-0.5,+0.5)
fs=1;                       %sampling frequency (Hz)
Nseg=1; L=N/Nseg;           %one full-length segment
[pxx1,f1]=pwelch(x,hamming(L),L/2,L,fs);
plot(f1,pxx1,'-b'); xlabel('Frequency (Hz)'); ylabel('Pxx'); grid on; hold on

The lines above create a signal, and compute the spectrum with one full-length segment. It is like your code, except I have used a Hamming window to taper the data. This is not necessary for a random sequence, but it is a good idea for spectrum estimation from real experimental data.

Nseg=16; L=N/Nseg;
[pxx16,f16]=pwelch(x,hamming(L),L/2,L,fs);
plot(f16,pxx16,'-r'); legend('Nseg=1','Nseg=16');

The lines above compute the spectrum by dividing the signal into 31 half-overlapped segments (31=2*Nseg-1). Each segment is 1/16 times as long as the original signal. The spectra are averaged to reduce the variance of the spectrum estimate. Next: repeat with Nseg=128 and Nseg=1024, to reduce the variance of the estimate even more.

Nseg=128; L=N/Nseg;
[pxx128,f128]=pwelch(x,hamming(L),L/2,L,fs);
Nseg=1024; L=N/Nseg;
[pxx1k,f1k]=pwelch(x,hamming(L),L/2,L,fs);

Make another plot, with the estimates with Nseg=16,128,1024. The plot shows that the variance decreases, and the spectral resolution decreases, as Nseg increases. Notice that the vertical axis range of the plot below is very different from the vertical axis range of the plot above. The red data plotted below is the same red data that was plotted above.

figure; plot(f16,pxx16,'-r',f128,pxx128,'-g.',f1k,pxx1k,'-b.'); 
legend('Nseg=16','Nseg=128','Nseg=1024'); grid on

Good luck.

[William Rose]

@Vlad Atanasiu,

The expected value for the one-sided spectrum of a random signal with zero mean is

where is the variance of the time-domain random signal. In this case, the variance of U([-0.5,+0.5) is 1/12, so the expected value for the spectrum is 1/6. The plots above are consistent with this prediction. If you use Gaussian noise, i.e.

x=randn(1,N);

then , and the expected value for the one-sided spectrum is 2*1=2. If you estimate the spectrum of such a signal, using the methods above, you will confirm this prediction.

[Vlad Atanasiu]

Thank you, William, for your detailed reply. I sum up below the answer to my question, focusing on the issue of the shape of the power spectrum of uniform random white noise.

1. A single uniform white noise signal of finite length will not have a flat power spectrum as expected. This is so because the amounts of samples from which the power at each frequency is computed are not equal. The observed non-uniform distribution is an artifact of the finite signal length.

2. To obtain the expected flat spectrum for uniform random white noise it is needed to either (a) average the spectra of segments of a single such signal, or (b) average the spectra of multiple such signals. The averaging ensures that there is an equal number of samples per frequency.

3. Due to the averaging process, the amplitude of an averaged power spectrum of uniform random white noise is not uniform random, but a normal distribution in the limit, as per the central limit theorem. Hence the Gaussian distribution observed in the diagrams.

Averaging the power spectra of multiple signals

Averaging the power spectra of segments of a single signal

% Investigating the shape and amplitude distribution of the power spectrum
% of a signal with uniform random amplitude values.

%% AVERAGING THE POWER SPECTRA OF MULTIPLE SIGNALS

clear all

% signal
k = 10^2; % number of signal realizations
n = 2^20; % 2^9 = 512, 2^13 = 8192 % signal length
t = 1:n; % time
y = unifrnd(-1,1,[n,k]); % uniform random amplitude

% spectral analysis
nyq = ceil(n/2); % Nyquist frequency
w = 1:nyq; % frequency
f = fft(y); % Fourier transform
f = f(1:nyq,:); % half the symmetrical spectrum
m = abs(f); % magnitude
s = m.^2; % power
a = (mean(s,2))'; % averaged power over all signals

% sorted values
ys = (sort(y(:,1),'descend'))';
as = sort(a,'descend');

% histogram
hn = 100;
hx = 1:hn;
ha = histcounts(a,hn);

% visualization
figure('Name','Analysis of noise spectra')

% signal
subplot(2,3,1)
plot(t(1:100),y(1:100,1),'b')
xlim('padded')
ylim('padded')
axis square
box on
title({'Detail of single signal realization',...
    'with uniform random amplitude',...
    ['and total length ',num2str(n)]})
xlabel('Time')
ylabel('Amplitude')

subplot(2,3,4)
plot(t,ys,'b')
xlim('padded')
ylim('padded')
axis square
box on
title({'Sorted amplitude values',...
    'of the entire signal'})
xlabel('Samples')
ylabel('Amplitude')

% averaged power spectrum
subplot(2,3,2)
plot(w,a,'b')
xlim('padded')
ylim('padded')
axis square
box on
title({'Averaged power spectra',['of ',num2str(k),' signals']})
xlabel('Frequency')
ylabel('Amplitude')

subplot(2,3,5)
plot(w,as,'b')
xlim('padded')
ylim('padded')
axis square
box on
title({'Sorted amplitude values',...
    'of the averaged power spectrum'})
xlabel('Samples')
ylabel('Amplitude')

% histogram of the power spectrum
subplot(2,3,3)
stairs(ha,hx,'b')
xlim('padded')
ylim('padded')
axis square
title({'Histogram','of averaged power spectrum'})
xlabel('Counts')
ylabel('Bins')



%% AVERAGING THE POWER SPECTRA OF SEGMENTS OF A SINGLE SIGNAL

clear all

% signal
n = 2^20; % 2^9 = 512, 2^13 = 8192 % signal length
t = 1:n; % time
y = unifrnd(-1,1,[n,1]); % uniform random amplitude

% spectral averaging using Welch’s power spectral density estimate
k = 100; % number of windows
nsc = floor(n/k);
nov = floor(nsc/2);
nff = max(256,2^nextpow2(nsc));
s = pwelch(y,hamming(nsc),nov,nff);
ns = size(s,1);
ws = 1:ns;

% sorted values
ys = (sort(y(:,1),'descend'))';
as = sort(s,'descend');

% histogram
hn = 100;
hx = 1:hn;
ha = histcounts(s,hn);

% visualization
figure('Name','Analysis of noise spectra')

% signal
subplot(2,3,1)
plot(t(1:100),y(1:100,1),'b')
xlim('padded')
ylim('padded')
axis square
box on
title({'Detail of single signal realization',...
    'with uniform random amplitude',...
    ['and total length ',num2str(n)]})
xlabel('Time')
ylabel('Amplitude')

subplot(2,3,4)
plot(t,ys,'b')
xlim('padded')
ylim('padded')
axis square
box on
title({'Sorted amplitude values',...
    'of the entire signal'})
xlabel('Samples')
ylabel('Amplitude')

% averaged power spectrum
subplot(2,3,2)
plot(ws,s,'b')
xlim('padded')
ylim('padded')
axis square
box on
title({'Averaged power spectrum',...
    ['of ',num2str(2*floor(k)),' signal segments']})
xlabel('Frequency')
ylabel('Amplitude')

subplot(2,3,5)
plot(ws,as,'b')
xlim('padded')
ylim('padded')
axis square
box on
title({'Sorted amplitude values',...
    'of the averaged power spectrum'})
xlabel('Samples')
ylabel('Amplitude')

% histogram of the power spectrum
subplot(2,3,3)
stairs(ha,hx,'b')
xlim('padded')
ylim('padded')
axis square
title({'Histogram','of averaged power spectrum'})
xlabel('Counts')
ylabel('Bins')

[William Rose]

You're welcome, @Vlad Atanasiu.