www.gusucode.com > signal 案例源码程序 matlab代码 > signal/ProbabilityOfDetectionAndExponentialDistributionExample.m

    %% Exponential Distribution
% Generate $N=1024$ samples of white noise with variance $\sigma=1$, given
% a sample rate of 1 Hz. Compute the power spectrum of the white noise.
% Choose the Lomb-Scargle normalization and specify an oversampling factor
% ${\tt ofac}=15$. Reset the random number generator for repeatable
% results.

% Copyright 2015 The MathWorks, Inc.


rng default

N = 1024;
t = (1:N)';
wgn = randn(N,1);

ofac = 15;
[pwgn,f] = plomb(wgn,t,[],ofac,'normalized');

%%
% Verify that the Lomb-Scargle power spectrum estimate of white noise has
% an exponential distribution with unit mean. Plot a histogram of the
% values of |pwgn| and a histogram of a set of exponentially distributed
% random numbers generated using the distribution function $f(z|1)=e^{-z}$.
% To normalize the histograms, recall that the total number of periodogram
% samples is $N\times{\tt ofac}/2$. Specify a bin width of 0.25. Overlay a
% plot of the theoretical distribution function.

dx = 0.25;
br = 0:dx:5;

Nf = N*ofac/2;

hpwgn = histcounts(pwgn,br)';

hRand = histcounts(-log(rand(Nf,1)),br)';

bend = br(1:end-1);

bar(bend,[hpwgn hRand]/Nf/dx,'histc')
hold on
plot(br,exp(-br))
legend('wgn','Empirical pdf','Theoretical pdf')
hold off

%%
% Embed in the noise a sinusoidal signal of frequency 0.1 Hz. Use a
% signal-to-noise ratio of $\xi=0.01$. Specify the sinusoid amplitude,
% $x_0$, using the relation $x_0=\sigma\sqrt{2\xi}$. Compute the power
% spectrum of the signal and plot its histogram alongside the empirical and
% theoretical distribution functions.

SNR = 0.01;
x0 = sqrt(2*SNR);
sigsmall = wgn + x0*sin(2*pi/10*t);

[psigsmall,f] = plomb(sigsmall,t,[],ofac,'normalized');

hpsigsmall = histcounts(psigsmall,br)';

bar(bend,[hpsigsmall hRand]/Nf/dx,'histc')
hold on
plot(br,exp(-br))
legend('sigsmall','Empirical pdf','Theoretical pdf')
hold off

%%
% Repeat the calculation using $\xi=1$. The distribution now differs
% noticeably.

SNR = 1;
x0 = sqrt(2*SNR);
siglarge = wgn + x0*sin(2*pi/10*t);

[psiglarge,f] = plomb(siglarge,t,[],ofac,'normalized');

hpsiglarge = histcounts(psiglarge,br)';

bar(bend,[hpsiglarge hRand]/Nf/dx,'histc')
hold on
plot(br,exp(-br))
legend('siglarge','Empirical pdf','Theoretical pdf')
hold off