modulate

function y=modulate(x,fs,f0)
%  y=modulate(x,fs,f0)
% y(t)=cos(2*pi*f0*t)*x(t)
    y=cos((2*pi*f0/fs)*[0:length(x)-1]).*x;