function[psd]=mpskfskpsd(T,fc,df,K,dt,M)

%T=Tb
%K=Number of Symbols
%dt=Number of sample points per symbol
%M=Number of runs
%fc=Center Frequency
%df=Frequency deviation of f1 and f2 to fc

t=0:dt*T:T-T*dt;
w=0;
for m=1:M,
    s1=zeros([1 K/dt]);
    s2=s1;
    s3=s1;
    for k=0:2:K-1,
        b=sign(randn([1 2]));
        q=(1+b)/2*[1 2]';
        Dt=k/dt+1:(k+2)/dt;
        s1(Dt)=cos([2*pi*(fc+b(1)*df)*(k*T+t) 2*pi*(fc+b(2)*df)*((k+1)*T+t)]);
        s2(Dt)=[b(1)*cos(2*pi*fc*(k*T+t)) b(2)*cos(2*pi*fc*((k+1)*T+t))];
        s3(Dt)=2^0.5*cos(2*pi*fc*([k*T+t (k+1)*T+t])+q*pi/2+pi/4);
        end    
    ps=abs(fftshift(fft([s1;s2;s3]')'*dt*T)).^2;
    w=w+ps/(K*T);
end

psd=w/M;
plot(-1/(2*dt*T):1/(K*T):1/(2*dt*T)-1/(K*T),10*log10(psd)),grid on