Covariance matrix $P$ and $Q$ for complex valued data The equations for the Kalman filtering are

  IM = H*X;
  IS = (R + H*P*H');
  K = P*H'/IS;
  X = X + K * (y-IM);
  P = P - K*IS*K';


The covariance matrix in my application is of dimension 3 by 3 size because the there are 3 state variables that need to be estimated.

For data in the real domain, I initialized the covariance matrices $Q$ (process covaraince) and $P$ (prediction covariace) in Matlab as

P= 5*(eye(3,3))
Q  = 0.001*(eye(3,3));

For the complex domain, I initialized $P$ and $Q$ as

    P_real = 5*(eye(d,d));
    P_imag = 5*(eye(d,d));
    Q_real = eye(d,d)* 10^(-2);
    Q_imag = eye(d,d)* 10^(-2);

    P  = P_real + sqrt(-1)*P_imag;

    Q  = Q_real + sqrt(-1)*Q_imag;


R is a scalar real valued number.  

Problem  : The default operation in Matlab is complex conjugate represented by the operator (.) '. I am using a toolbox which applies the above expression for Kalman filtering and there is no term conj() present especially when evaluating the expression IS , K, P. However, the toolbox never uses any (.).' or (.).* operators as well. When in complex domain, would the equations change to

IM = H*X  % no conjugate

IS = (R + conj(H)*P*ctranspose(H))

K = P*ctranspose(H)/IS;

X = X + K *(y - IM)  % here * denotes the multiplication operation

P = P - conj(K)*IS*ctraspose(K)

where H, X, y, P are all complex valued.

Are these okay?
In general, the equations from your first code block work equally for real and complex values. 

Note that the ()' operation is the Hermitian Conjugate (i.e. transpose + conjugate). If you are in the real domain, it only becomes transpose (because conjugation on reals doesn't change them). Hence, you would not need to change anything in the code. 

Regarding your covariance matrices: Note that a covariance matrix is always positive (semi-)definite and Hermitian symmetric, i.e. $P^H=P$. This ensures that each single variance on the diagonal is a real value. Also, how would you imagine a e.g. Gaussian with complex variance?

However, there is one thing which might be important for you here: When considering white Gaussian noise in the real domain, you generate random values according to

noise_real = sqrt(sigma2) * randn(N,M);


where sigma2 is the variance of this random variable. If you are in the complex domain, you have to generate them via

noise_complex = sqrt(sigma2/2) * (randn(N,M) + 1j*randn(N,M))


i.e. the variance becomes halfed for real and imaginary domain. Then, when you calculate var(noise_complex) you will get exactly sigma2.

I see from your post that your P and Q denote uncorrelated complex variables. In this case, if they are supposed to have a variance of sigma2 on each element, this is the corresponding code:

P = sigma2*eye(3);