fs=10^5; N=10^5 ;
t=1/fs:1/fs:N/fs;
fm=2; fc=10^3;
m= sin(2*pi*fm*t);
c= sin(2*pi*fc*t);
v= m.*c;

SNRdB =10; SNR = 10^(SNRdB/10);
vn= var (v)/SNR;
n= sqrt(vn)*randn(1,N);
v= v+n;

r=zeros(1,N); n=fs/fc;
for k=1:fc
    mr ((k-1)*n+1:k*n)=2*v((k-1)*n+1:k*n)*c((k-1)*n+1:k*n)'/n;
end
figure(1)
subplot(2,1,1); plot(t,m);
subplot(2,1,2); plot(t,mr);axis([0 1 -1 1])