function [I,Q]=mod_16QAM(x)
  
% [I,Q]=modulate(x,modulation)
%
% 16QAM modulation
% x : input bit stream
% modulation : 'bpsk','qpsk','8psk','16qam', or '64qam'
% [I,Q] : complex symbols
% x,I,Q are column vectors

i=sqrt(-1);

if rem(length(x),4)
    error('Length of input vector must be divisible with 4 for 16qam');
end
b0=x(1:4:end); b1=x(2:4:end); b2=x(3:4:end); b3=x(4:4:end);
scale=3/sqrt(10);
y=((-1).^b0.*(1-2/3*b1) + i*(-1).^b2.*(1-2/3*b3))*scale;

I=real(y);
Q=imag(y);