Reply by ryujin_ssdt July 24, 20072007-07-24
>On Jul 24, 4:52 am, "ryujin_ssdt" <hsan...@gmail.com> wrote: >> I read that the spectrum of a baseband BPSK signal, if the symbols are >> equiproblables and uncorrelated, is the same as the spectrum of the >> shaping pulse used. >> >> I tried a small simulation to visualize this but I think I am missing >> something. I use a Raised Cosine filter (following matlabs own
examples)
>> for pulse shaping a random generated BPSK signal. >> >> I plot the theorical spectrum of the raised cosine (taken from Proakis >> book) and the simulated spectrum of the random BPSK signal using
matlab's
>> periodogram method and by direct FFT method. >> >> So far the simulated graphs have nothing similar to the theory.
Increasing
>> the number of samples (more randomness) seems to slowly approximate to
the
>> theory but still even for 3e4 points the simulated results are far
from
>> the theorical one. >> >> Any thoughts on these results? surely I am doing something wrong or >> missing an important concept. >> >> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> % BPSK PSD estimation test >> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> clear; clc; >> >> M = 2; % Alphabet Size >> k = log2(M); % Number of Bits per symbol >> fs = 1; % Sampling frequency >> T = 1/fs; % Sampling period >> t = [0:T:100]; % time vector >> nsamp = 5; % Oversampling rate >> fftlength = 2^nextpow2(length(t)); % FFT length >> ebno = 5; % Eb/No >> >> % Initialize the root-raised cosine (Pulse Shaping Filter) >> filtorder = 20; % Filter order >> delay = filtorder/(nsamp*2); % Group delay (# of input samples) >> rolloff = 0.5; % Rolloff factor of filter >> >> % Create a square root raised cosine filter. >> rrcfilter = rcosine(1,nsamp,'fir/sqrt',rolloff,delay); >> >> % Signal Source >> x = randint(length(t),1); % Random binary data >> >> % Now transmit the signal throught a channel with >> % different EbNo values. >> >> % Modulate the signal >> ytx = pskmod(x,M); >> >> % Pulse Shape and oversample the signal >> ytx = rcosflt(ytx,1,nsamp,'filter',rrcfilter); >> >> % add some AWGN noise >> snr = ebno + 10*log10(k) - 10*log10(nsamp); >> ynoisy = awgn(ytx, snr, 'measured'); >> >> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> %% Theorical BPSK PSD (Raised Cosine) >> %% from Proakis - Digital Communications 4th edition pg 547 >> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> f = [0:1/fftlength:fs/2]; >> Px = T/2 * (1 + cos((pi*T/rolloff)*(abs(f)-((1 - rolloff)/(2*T))))); >> Px(f < (1-rolloff)/(2*T)) = T; >> Px(f > (1+rolloff)/(2*T)) = 0; >> >> figure(1); >> grid on; >> subplot(311); >> plot(f,Px) >> title('Theorical BPSK PSD'); >> xlabel('Hz'); >> >> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> %% PSD estimation via Matlab periodogram method >> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> W = blackman(length(ynoisy)); >> Px = periodogram(ynoisy, W, fftlength, fs); >> subplot(312); >> f = [0:1/fftlength:fs/2]; >> plot(f,Px(1:length(f)), 'r'); >> title('Periodogram Method'); >> xlabel('Hz'); >> >> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> %% PSD estimation via direct FFT >> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> X = fft(ynoisy.*W,fftlength); >> f = [0:1/fftlength:fs/2]; >> Px=X.*conj(X); >> >> subplot(313); >> plot(f,Px(1:length(f)), 'r'); >> title('Estimated BPSK PSD (FFT method)'); >> xlabel('Hz'); > >All you're really testing here is whether or not applying a RRC filter >to a zero-stuffed series of +/- impulses will generate a PSD of the >theoretical shape, not really anything specific to BPSK. I can see a >couple reasons why you're not getting exactly what you expect: >
>- Eb/No is set to 5 dB in your script, then later when calculating the >SNR of the noisy signal, you subtract 10*log10(nsamp) from it, which >gives you a negative SNR. Thus, the result looks more like junk. Try >increasing Eb/No to get nicer plots.
Indeed increasing Eb/No or even removing the noise altogether gives better results but still far away from theory.
>- Since your data is oversampled 5x, the resulting spectrum that you >get from the FFT or periodogram will be frequency-compressed into one- >fifth of the full range. If you account for the oversampling in your >theoretical calculation, the curves will look a lot more similar.
I think it will take a little longer for me to digest this. I don't understand yet the implications of that oversampling (just copied from examples) and if possible I would like to remove it (nsamp = 1) but if I do set nsamp = 1, then the rcosine method complains. Any document references on why this oversampling is needed would be appreciated.
>- You might also try using dB scales for the plots; they won't look as >jagged as the linear-scaled ones you have now.
This one too improved the resulting graphs. At least I can see something familiar to a sinc function but still far from theory. I was hoping to see something like the BER graphs were the theoretical and simulated curves almost coincide but guess that wont be the case here. To convert to dB I simply did 10*log10(Px). Is this the correct way of doing this??
> >Jason > >
Below modified source %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % BPSK PSD estimation test %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear; clc; M = 2; % Alphabet Size k = log2(M); % Number of Bits per symbol fs = 1; % Sampling frequency T = 1/fs; % Sampling period t = [0:T:100]; % time vector nsamp = 2; % Oversampling rate fftlength = 2^nextpow2(length(t)); % FFT length ebno = 25; % Eb/No % Initialize the root-raised cosine (Pulse Shaping Filter) filtorder = 20; % Filter order delay = filtorder/(nsamp*2); % Group delay (# of input samples) rolloff = 0.5; % Rolloff factor of filter % Create a square root raised cosine filter. rrcfilter = rcosine(1,nsamp,'fir/sqrt',rolloff,delay); % Signal Source x = randint(length(t),1); % Random binary data % Now transmit the signal throught a channel with % different EbNo values. % Modulate the signal ytx = pskmod(x,M); % Pulse Shape and oversample the signal ytx = rcosflt(ytx,1,nsamp,'filter',rrcfilter); % add some AWGN noise snr = ebno + 10*log10(k) - 10*log10(nsamp); ynoisy = awgn(ytx, snr, 'measured'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Theorical BPSK PSD (Raised Cosine) %% from Proakis - Digital Communications 4th edition pg 547 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% f = [0:fftlength/2]*fs/fftlength; Px = T/2 * (1 + cos((pi*T/rolloff)*(abs(f)-((1 - rolloff)/(2*T))))); Px(f < (1-rolloff)/(2*T)) = T; Px(f > (1+rolloff)/(2*T)) = 0; figure(1); grid on; subplot(311); plot(f,10*log10(Px)) title('Theorical BPSK PSD'); xlabel('Hz'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PSD estimation via Matlab periodogram method %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% W = blackman(length(ynoisy)); Px = periodogram(ynoisy, W, fftlength, fs); subplot(312); f = [0:fftlength/2]*fs/fftlength; plot(f,10*log10(Px(1:length(f))), 'r'); title('Periodogram Method'); xlabel('Hz'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% PSD estimation via Periodogram %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% X = fft(ynoisy.*W,fftlength); f = [0:fftlength/2]*fs/fftlength; Px=X.*conj(X)./length(ynoisy); subplot(313); plot(f,10*log10(Px(1:length(f))), 'r'); title('Estimated BPSK PSD (FFT method)'); xlabel('Hz');
Reply by ryujin_ssdt July 24, 20072007-07-24
>On Jul 24, 4:52 am, "ryujin_ssdt" <hsan...@gmail.com> wrote: >> I read that the spectrum of a baseband BPSK signal, if the symbols are >> equiproblables and uncorrelated, is the same as the spectrum of the >> shaping pulse used. >> >> I tried a small simulation to visualize this but I think I am missing >> something. I use a Raised Cosine filter (following matlabs own
examples)
>> for pulse shaping a random generated BPSK signal. >> >> I plot the theorical spectrum of the raised cosine (taken from Proakis >> book) and the simulated spectrum of the random BPSK signal using
matlab's
>> periodogram method and by direct FFT method. >> >> So far the simulated graphs have nothing similar to the theory.
Increasing
>> the number of samples (more randomness) seems to slowly approximate to
the
>> theory but still even for 3e4 points the simulated results are far
from
>> the theorical one. >> >> Any thoughts on these results? surely I am doing something wrong or >> missing an important concept. >> >> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> % BPSK PSD estimation test >> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> clear; clc; >> >> M = 2; % Alphabet Size >> k = log2(M); % Number of Bits per symbol >> fs = 1; % Sampling frequency >> T = 1/fs; % Sampling period >> t = [0:T:100]; % time vector >> nsamp = 5; % Oversampling rate >> fftlength = 2^nextpow2(length(t)); % FFT length >> ebno = 5; % Eb/No >> >> % Initialize the root-raised cosine (Pulse Shaping Filter) >> filtorder = 20; % Filter order >> delay = filtorder/(nsamp*2); % Group delay (# of input samples) >> rolloff = 0.5; % Rolloff factor of filter >> >> % Create a square root raised cosine filter. >> rrcfilter = rcosine(1,nsamp,'fir/sqrt',rolloff,delay); >> >> % Signal Source >> x = randint(length(t),1); % Random binary data >> >> % Now transmit the signal throught a channel with >> % different EbNo values. >> >> % Modulate the signal >> ytx = pskmod(x,M); >> >> % Pulse Shape and oversample the signal >> ytx = rcosflt(ytx,1,nsamp,'filter',rrcfilter); >> >> % add some AWGN noise >> snr = ebno + 10*log10(k) - 10*log10(nsamp); >> ynoisy = awgn(ytx, snr, 'measured'); >> >> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> %% Theorical BPSK PSD (Raised Cosine) >> %% from Proakis - Digital Communications 4th edition pg 547 >> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> f = [0:1/fftlength:fs/2]; >> Px = T/2 * (1 + cos((pi*T/rolloff)*(abs(f)-((1 - rolloff)/(2*T))))); >> Px(f < (1-rolloff)/(2*T)) = T; >> Px(f > (1+rolloff)/(2*T)) = 0; >> >> figure(1); >> grid on; >> subplot(311); >> plot(f,Px) >> title('Theorical BPSK PSD'); >> xlabel('Hz'); >> >> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> %% PSD estimation via Matlab periodogram method >> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> W = blackman(length(ynoisy)); >> Px = periodogram(ynoisy, W, fftlength, fs); >> subplot(312); >> f = [0:1/fftlength:fs/2]; >> plot(f,Px(1:length(f)), 'r'); >> title('Periodogram Method'); >> xlabel('Hz'); >> >> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> %% PSD estimation via direct FFT >> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> X = fft(ynoisy.*W,fftlength); >> f = [0:1/fftlength:fs/2]; >> Px=X.*conj(X); >> >> subplot(313); >> plot(f,Px(1:length(f)), 'r'); >> title('Estimated BPSK PSD (FFT method)'); >> xlabel('Hz'); > > > >The RRC is normally applied *before* modulation. > >John > >
Thanks John for your fast comment and you are correct in that RRC is done before modulation but I think when you say "modulation" that you mean "Bandpass modulation". If you see I am working only with Baseband BPSK (i.e. no carriers). You can verify the Matlab example (edit commdoc_rrc) from the Communication Toolbox that I used as reference to implement my script. thanks Horacio
Reply by July 24, 20072007-07-24
On Jul 24, 4:52 am, "ryujin_ssdt" <hsan...@gmail.com> wrote:
> I read that the spectrum of a baseband BPSK signal, if the symbols are > equiproblables and uncorrelated, is the same as the spectrum of the > shaping pulse used. > > I tried a small simulation to visualize this but I think I am missing > something. I use a Raised Cosine filter (following matlabs own examples) > for pulse shaping a random generated BPSK signal. > > I plot the theorical spectrum of the raised cosine (taken from Proakis > book) and the simulated spectrum of the random BPSK signal using matlab's > periodogram method and by direct FFT method. > > So far the simulated graphs have nothing similar to the theory. Increasing > the number of samples (more randomness) seems to slowly approximate to the > theory but still even for 3e4 points the simulated results are far from > the theorical one. > > Any thoughts on these results? surely I am doing something wrong or > missing an important concept. > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > % BPSK PSD estimation test > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > clear; clc; > > M = 2; % Alphabet Size > k = log2(M); % Number of Bits per symbol > fs = 1; % Sampling frequency > T = 1/fs; % Sampling period > t = [0:T:100]; % time vector > nsamp = 5; % Oversampling rate > fftlength = 2^nextpow2(length(t)); % FFT length > ebno = 5; % Eb/No > > % Initialize the root-raised cosine (Pulse Shaping Filter) > filtorder = 20; % Filter order > delay = filtorder/(nsamp*2); % Group delay (# of input samples) > rolloff = 0.5; % Rolloff factor of filter > > % Create a square root raised cosine filter. > rrcfilter = rcosine(1,nsamp,'fir/sqrt',rolloff,delay); > > % Signal Source > x = randint(length(t),1); % Random binary data > > % Now transmit the signal throught a channel with > % different EbNo values. > > % Modulate the signal > ytx = pskmod(x,M); > > % Pulse Shape and oversample the signal > ytx = rcosflt(ytx,1,nsamp,'filter',rrcfilter); > > % add some AWGN noise > snr = ebno + 10*log10(k) - 10*log10(nsamp); > ynoisy = awgn(ytx, snr, 'measured'); > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > %% Theorical BPSK PSD (Raised Cosine) > %% from Proakis - Digital Communications 4th edition pg 547 > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > f = [0:1/fftlength:fs/2]; > Px = T/2 * (1 + cos((pi*T/rolloff)*(abs(f)-((1 - rolloff)/(2*T))))); > Px(f < (1-rolloff)/(2*T)) = T; > Px(f > (1+rolloff)/(2*T)) = 0; > > figure(1); > grid on; > subplot(311); > plot(f,Px) > title('Theorical BPSK PSD'); > xlabel('Hz'); > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > %% PSD estimation via Matlab periodogram method > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > W = blackman(length(ynoisy)); > Px = periodogram(ynoisy, W, fftlength, fs); > subplot(312); > f = [0:1/fftlength:fs/2]; > plot(f,Px(1:length(f)), 'r'); > title('Periodogram Method'); > xlabel('Hz'); > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > %% PSD estimation via direct FFT > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > X = fft(ynoisy.*W,fftlength); > f = [0:1/fftlength:fs/2]; > Px=X.*conj(X); > > subplot(313); > plot(f,Px(1:length(f)), 'r'); > title('Estimated BPSK PSD (FFT method)'); > xlabel('Hz');
All you're really testing here is whether or not applying a RRC filter to a zero-stuffed series of +/- impulses will generate a PSD of the theoretical shape, not really anything specific to BPSK. I can see a couple reasons why you're not getting exactly what you expect: - Eb/No is set to 5 dB in your script, then later when calculating the SNR of the noisy signal, you subtract 10*log10(nsamp) from it, which gives you a negative SNR. Thus, the result looks more like junk. Try increasing Eb/No to get nicer plots. - Since your data is oversampled 5x, the resulting spectrum that you get from the FFT or periodogram will be frequency-compressed into one- fifth of the full range. If you account for the oversampling in your theoretical calculation, the curves will look a lot more similar. - You might also try using dB scales for the plots; they won't look as jagged as the linear-scaled ones you have now. Jason
Reply by John July 24, 20072007-07-24
On Jul 24, 4:52 am, "ryujin_ssdt" <hsan...@gmail.com> wrote:
> I read that the spectrum of a baseband BPSK signal, if the symbols are > equiproblables and uncorrelated, is the same as the spectrum of the > shaping pulse used. > > I tried a small simulation to visualize this but I think I am missing > something. I use a Raised Cosine filter (following matlabs own examples) > for pulse shaping a random generated BPSK signal. > > I plot the theorical spectrum of the raised cosine (taken from Proakis > book) and the simulated spectrum of the random BPSK signal using matlab's > periodogram method and by direct FFT method. > > So far the simulated graphs have nothing similar to the theory. Increasing > the number of samples (more randomness) seems to slowly approximate to the > theory but still even for 3e4 points the simulated results are far from > the theorical one. > > Any thoughts on these results? surely I am doing something wrong or > missing an important concept. > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > % BPSK PSD estimation test > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > clear; clc; > > M = 2; % Alphabet Size > k = log2(M); % Number of Bits per symbol > fs = 1; % Sampling frequency > T = 1/fs; % Sampling period > t = [0:T:100]; % time vector > nsamp = 5; % Oversampling rate > fftlength = 2^nextpow2(length(t)); % FFT length > ebno = 5; % Eb/No > > % Initialize the root-raised cosine (Pulse Shaping Filter) > filtorder = 20; % Filter order > delay = filtorder/(nsamp*2); % Group delay (# of input samples) > rolloff = 0.5; % Rolloff factor of filter > > % Create a square root raised cosine filter. > rrcfilter = rcosine(1,nsamp,'fir/sqrt',rolloff,delay); > > % Signal Source > x = randint(length(t),1); % Random binary data > > % Now transmit the signal throught a channel with > % different EbNo values. > > % Modulate the signal > ytx = pskmod(x,M); > > % Pulse Shape and oversample the signal > ytx = rcosflt(ytx,1,nsamp,'filter',rrcfilter); > > % add some AWGN noise > snr = ebno + 10*log10(k) - 10*log10(nsamp); > ynoisy = awgn(ytx, snr, 'measured'); > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > %% Theorical BPSK PSD (Raised Cosine) > %% from Proakis - Digital Communications 4th edition pg 547 > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > f = [0:1/fftlength:fs/2]; > Px = T/2 * (1 + cos((pi*T/rolloff)*(abs(f)-((1 - rolloff)/(2*T))))); > Px(f < (1-rolloff)/(2*T)) = T; > Px(f > (1+rolloff)/(2*T)) = 0; > > figure(1); > grid on; > subplot(311); > plot(f,Px) > title('Theorical BPSK PSD'); > xlabel('Hz'); > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > %% PSD estimation via Matlab periodogram method > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > W = blackman(length(ynoisy)); > Px = periodogram(ynoisy, W, fftlength, fs); > subplot(312); > f = [0:1/fftlength:fs/2]; > plot(f,Px(1:length(f)), 'r'); > title('Periodogram Method'); > xlabel('Hz'); > > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > %% PSD estimation via direct FFT > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% > X = fft(ynoisy.*W,fftlength); > f = [0:1/fftlength:fs/2]; > Px=X.*conj(X); > > subplot(313); > plot(f,Px(1:length(f)), 'r'); > title('Estimated BPSK PSD (FFT method)'); > xlabel('Hz');
The RRC is normally applied *before* modulation. John
Reply by ryujin_ssdt July 24, 20072007-07-24
I read that the spectrum of a baseband BPSK signal, if the symbols are
equiproblables and uncorrelated, is the same as the spectrum of the
shaping pulse used.

I tried a small simulation to visualize this but I think I am missing
something. I use a Raised Cosine filter (following matlabs own examples)
for pulse shaping a random generated BPSK signal.

I plot the theorical spectrum of the raised cosine (taken from Proakis
book) and the simulated spectrum of the random BPSK signal using matlab's
periodogram method and by direct FFT method.

So far the simulated graphs have nothing similar to the theory. Increasing
the number of samples (more randomness) seems to slowly approximate to the
theory but still even for 3e4 points the simulated results are far from
the theorical one.

Any thoughts on these results? surely I am doing something wrong or
missing an important concept.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BPSK PSD estimation test
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clc;

M = 2;                              % Alphabet Size
k = log2(M);                        % Number of Bits per symbol
fs = 1;                             % Sampling frequency
T = 1/fs;                           % Sampling period
t = [0:T:100];                      % time vector
nsamp = 5;                          % Oversampling rate
fftlength = 2^nextpow2(length(t));  % FFT length
ebno = 5;                           % Eb/No

% Initialize the root-raised cosine (Pulse Shaping Filter)
filtorder = 20;                 % Filter order
delay = filtorder/(nsamp*2);    % Group delay (# of input samples)
rolloff = 0.5;                  % Rolloff factor of filter

% Create a square root raised cosine filter.
rrcfilter = rcosine(1,nsamp,'fir/sqrt',rolloff,delay);

% Signal Source
x = randint(length(t),1);   % Random binary data

% Now transmit the signal throught a channel with
% different EbNo values.

% Modulate the signal
ytx = pskmod(x,M);

% Pulse Shape and oversample the signal
ytx = rcosflt(ytx,1,nsamp,'filter',rrcfilter);

% add some AWGN noise
snr = ebno + 10*log10(k) - 10*log10(nsamp);
ynoisy = awgn(ytx, snr, 'measured');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Theorical BPSK PSD (Raised Cosine)
%% from Proakis - Digital Communications 4th edition pg 547
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f = [0:1/fftlength:fs/2];
Px = T/2 * (1 + cos((pi*T/rolloff)*(abs(f)-((1 - rolloff)/(2*T)))));
Px(f < (1-rolloff)/(2*T)) = T;
Px(f > (1+rolloff)/(2*T)) = 0;

figure(1);
grid on;
subplot(311);
plot(f,Px)
title('Theorical BPSK PSD');
xlabel('Hz');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PSD estimation via Matlab periodogram method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
W = blackman(length(ynoisy));
Px = periodogram(ynoisy, W, fftlength, fs);
subplot(312);
f = [0:1/fftlength:fs/2];
plot(f,Px(1:length(f)), 'r');
title('Periodogram Method');
xlabel('Hz');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PSD estimation via direct FFT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X = fft(ynoisy.*W,fftlength);
f = [0:1/fftlength:fs/2];
Px=X.*conj(X);

subplot(313);
plot(f,Px(1:length(f)), 'r');
title('Estimated BPSK PSD (FFT method)');
xlabel('Hz');