Frequency Translation in MATLAB
You are partly right. However It's not the filter coefficinets $a$,$b$, but the filter impulse response $h[n]$ which you have to multiply with the exponential term $e^{j w_0 n}$ to shift its frequency response to a desired center of frequency $w_0$. Note that for FIR filters the coefficients are such that $a=1$ and $b[n]=h[n]$; i.e., the impulse response equals the coefficients $b[k]$, but that's not so with the IIR filters.

Assume that you have an FIR filter with coefficients $a=1$ and $b[k]$ so that its impulse response is: $$h[n] = b[n]$$ And let the frequency response of the filter be $H(e^{j \omega})$, then the frequency response of the new filter $h_+[n] = e^{j w_0 n} h[n]$ will be
$$ h_+[n] = e^{j w_0 n} h[n] \implies H_+(e^{j \omega}) = H(e^{j (\omega - \omega_0)}) $$

and similarly the frequency response of the new filter $h_-[n] = e^{-j w_0 n} h[n]$ will be
$$ h_-[n] = e^{-j \omega_0 n} h[n] \implies H_-(e^{j \omega}) = H_(e^{j( \omega + \omega_0)}) $$

Note if you want a real impulse response filter, say $h_0[n]$, then you should add those two sub shifted filters:
$$ h_0[n] = h_+[n] + h_-[n] \implies H_0(e^{j \omega}) = H_+(e^{j \omega}) + H_-(e^{j \omega}) = H(e^{j (\omega - \omega_0)}) + H(e^{j (\omega + \omega_0)}) $$

You can achieve the same result by the following multiplication then;
$$ h_0[n] = 2\cos(\omega_0 n) h[n] \implies H_0(e^{j \omega}) = H(e^{j (\omega - \omega_0)}) + H(e^{j (\omega + \omega_0)}) $$

MATLAB implmenetation would be as follows:

h = fir1(64, 0.1 );     % prototype lowpass filter
L = length(h);

w0 = 0.5*pi;            % bandpass filter new frequency
n = [0:L-1];            % discrete-time index vector 

h0 = 2*cos(w0*n).*h;    % new bandpass filter

% Display prototype lowpass FIR filter
figure,subplot(2,1,1)
stem(n,h);title('Prototype Lowpass filter h[n]');
subplot(2,1,2)
plot(linspace(-1,1,1024), abs(fftshift(fft(h,1024))));
title('DTFT magnitude |H(e^{j \omega})| of lowpass filter h[n]');

% Display new bandpass FIR filter
figure,subplot(2,1,1)
stem(n,h0);title('Prototype Lowpass filter h_0[n]');
subplot(2,1,2)
plot(linspace(-1,1,1024), abs(fftshift(fft(h0,1024))));
title('DTFT magnitude |H_0(e^{j \omega})| of bandpass filter h_0[n]');