%% Median Frequency of Sinusoids
% Generate 1024 samples of a 100.123 kHz sinusoid sampled at 1024 kHz. Add white Gaussian noise such that the signal-to-noise ratio is 40 dB. Reset the random number generator for reproducible results.

nSamp = 1024;
Fs = 1024e3;
SNR = 40;
rng default

t = (0:nSamp-1)'/Fs;

x = sin(2*pi*t*100.123e3);
x = x + randn(size(x))*std(x)/db2mag(SNR);

%
% Use the periodogram function to compute the power spectral density (PSD) of the signal. Specify a Kaiser window with the same length as the signal and a shape factor of 38. Estimate the median frequency of the signal and annotate it on a plot of the PSD.

[Pxx,f] = periodogram(x,kaiser(nSamp,38),[],Fs);

medfreq(Pxx,f);

%
% Generate another sinusoid, this one with a frequency of 257.321 kHz and an amplitude that is twice that of the first sinusoid. Add white noise.

x2 = 2*sin(2*pi*t*257.321e3);
x2 = x2 + randn(size(x2))*std(x2)/db2mag(SNR);

%
% Concatenate the sinusoids to produce a two-channel signal. Estimate the PSD of each channel and use the result to determine the median frequency.

[Pyy,f] = periodogram([x x2],kaiser(nSamp,38),[],Fs);

y = medfreq(Pyy,f)

%
% Annotate the median frequencies of the two channels on a plot of the PSDs.

medfreq(Pyy,f);

%
% Add the two channels to form a new signal. Estimate the PSD and annotate the median frequency.

[Pzz,f] = periodogram(x+x2,kaiser(nSamp,38),[],Fs);

medfreq(Pzz,f);

% Copyright 2015 The MathWorks, Inc.