I've said, in another comment, convolution using the function conv "i.e in MATLAB" and convolution using the Toeplitx matrix must give the same results. That's ok.
Now, according to your code, the received signal $r$ is the results of convolution between channel $h$ and emitted signal $x$, which means $r = h*x + n$
* indicate to the convolution (which is circular convolution in case of using $CP$ with OFDM or with any other system).
So, in that conventional known case, we are using SISO system, where 1 antenna is used as transmitter and 1 as receiver, The length of our parameters should be:
$x = N$ x 1 ; $h = L$ x 1;
$N$ is the length of our signal, in your case $Q$
$L$ the length of channel or we sometimes call it IR, in your case $M$
Till here, it's clear, it's the normal process which can be read anywhere. Now suppose that you are using $SIMO$ system with $P$ antennas at reciever instead of 1 (by the way fractional sampling is also equivalent to SIMO system). in that case you suppose to have $P$ copies for your signal instead of one. It's like you are doing $P$ times convolutions for your emitted signal with $P$ different channels. let's say $P$ = 4; means you have 4 receiver's antennas equivalent in our case into 4 different channels.
As mentioned, you suppose to have 4 copies for your signal, and 4 different channels compared with conventional case, so let's say parameter $H$ is toeplitz matrix which will represent those four channels in SIMO system. So we will have
$R = HX + N$
$R$ is received signal in SIMO, $H$ is the toeplitx matrix, $X$ emitted signal in SIMO system too and also $N$ represent the noise.
Now, your question how to build $H$ and what's dimension of $H$ when using SIMO system. (by they way, your code is correct)
The emitted signal $X$, represents 4 copies of $x$, so its dimension is $(N +L)$ x 1;
The toeplitz matrix $H$ has a dimension $NP$ x $(N+L)$ where $P$ = 4 in your case. And the noise $N$ should be $NP$ x 1. (N as noise is different about N of emitted signal in this case).
So, Now, you can conclude the dimension of received signal $r$ and $R$ in every case easily.
and regarding your question, why using $x1$, it's to have dimension of emitted signal $X$ of $(N +L)$ x 1 instead of dimension $N$ in conventional case