function bpskd(g,f)
if nargin > 2
error('Too many input arguments');
elseif nargin==1
    f=1;
end

if f<1;
    error('Frequency must be bigger than 1');
end

t=0:2*pi/99:2*pi;
cp=[];sp=[];
mod=[];mod1=[];bit=[];

for n=1:length(g);
    if g(n)==0; 
        die=-ones(1,100);   %Modulante
        se=zeros(1,100);    %Se馻l
    else g(n)==1;
        die=ones(1,100);    %Modulante
        se=ones(1,100);     %Se馻l
    end
    c=sin(f*t);
    cp=[cp die];    
    mod=[mod c];    
    bit=[bit se];
end

bpsk=cp.*mod;
subplot(4,1,1);plot(bit,'LineWidth',1.5);grid on;
title('Binary Signal');
axis([0 100*length(g) -2.5 2.5]);


subplot(4,1,2);plot(bpsk,'LineWidth',1.5);grid on;
title('Waveform of BPSK modulation in Time domain');
axis([0 100*length(g) -2.5 2.5]);

subplot(4,1,3);
SNR=20;
x=awgn(bpsk,SNR);
plot(x,'LineWidth',1.5);grid on;
title('Waveform of BPSK modulation with AWGN in Time domain');
axis([0 100*length(g) -2.5 2.5]);

subplot(4,1,4);
N=1024;
y=fft(x,N);%進行fft變換
mag=abs(y);%求幅值
f1=(0:length(y)-1)'*f/length(y);%進行對應的頻率轉換
plot(f1,mag);%做頻譜圖
xlabel('頻率(Hz)');
ylabel('幅值');
title('Waveform of BPSK modulation with AWGN in Frequency domain');
axis([0 f 0 300]);
%Py =2*(y.*conj(y))/N;
%計算功率譜密度Py
%figure;
%plot(f,Py);
%xlabel('頻率(Hz)');
%ylabel('功率譜密度');
%title('功率譜密度');
%grid;
%plot(y,'LineWidth',1.5);grid on;