Unexpected power spectrum of a signal with uniform random amplitude values

15 visualizaciones (últimos 30 días)
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
% 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')

Respuesta aceptada

William Rose
William Rose el 10 de En. de 2022
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.
  4 comentarios
Vlad Atanasiu
Vlad Atanasiu el 11 de En. de 2022
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
random_uniform_spectrum_analysis_averaging_multiple_signals
Averaging the power spectra of segments of a single signal
random_uniform_spectrum_analysis_averaging_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')

Iniciar sesión para comentar.

Más respuestas (0)

Productos


Versión

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by