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