fdesign.bandpass - Matching filters in Frequency and Time Domains

I am trying to use fdesign.bandpass for the first time.

I have a filter I have created manually. I wish to use fdesign to create a filter that will match the output of my manual filter in BOTH the frequency and time domain.

You will see that when I plot the below data post filtering, it clearly is not the same in the time domain. This is because a pair of filters with the same freq-domain response do not necessarily have the same time-domain response because they could be (non-linearly) phase-shifted relative to each other.

My ideal response to this question, would be if someone provide me with a modified version of the below code that will enforce time AND frequency domain outputs to be the same? Thank you!

function matlabQuestion_03()
%%http://www.mathworks.co.uk/help/dsp/ref/fdesign.bandpass.html
clc;
dbstop if error;
close all;

%%Generate Some Data
N = 5000;
data = cumsum(randn(N,1));

%%Create some Exponential filter coefficents
t_slow = 252;
alpha = 2 / (t_slow+1);
b = repmat(1-alpha,1 ,N).^(1:N); 
b = b ./ sum(b); 
a = 1;
t_fast = 100;
alpha = 2 / (t_fast+1);
c = repmat(1-alpha,1 ,N).^(1:N); 
c = c ./ sum(c); 
d = c-b;  %the difference of two exponential filters. simple.
maDiff = filter(d, a, data);

%%Plot the filters frequency response
%see this is bandpass-like with a long tail.
[h3 w3]= freqz(d,1, N); 
w3 = w3./pi;

figure;
subplot(211);
plot(w3, exp(abs(h3)),'g', 'linewidth', 3); % I am willing to chop off the long tail by specifiying an upper stop band.
title('I wish to recreate this filter using fdesign.bandpass');
xlim([0 0.3]);
subplot(212);
plot(log(w3), exp(abs(h3)));     

% get the "mid" frequency of the filter. How can i get this value from knowing t_slow and t_fast?
[~,idx] = max(exp(abs(h3)));
centerFrequency = w3(idx);

%%Now I would like to recreate maDiff using fdesign.bandpass
%"A bandpass filter can be formed by a parallel configuration of two low pass filters". p.624, Mandal and Asif (book). 

%place a spread around the mid frequency
Fc1 = 0.77^2 * centerFrequency; %Fc1 = lower_passBand
Fc2 = 1.23^2 * centerFrequency; %Fc2 = lower_passBand

Ast1 = 20;                   % attenuation in the first stop band in decibels (the default units). 
Ap   = 0.5;                  % amount of ripple allowed in the pass band.
Ast2 = 20;                   % attenuation in the second stop band in decibels (the default units). 
nOrder = 100;                %filter order for FIR filter
f = fdesign.bandpass('N,Fc1,Fc2,Ast1,Ap,Ast2',nOrder,Fc1,Fc2,Ast1,Ap,Ast2);
Hf = design(f);  
output = filter(Hf,data);   

%%Plot the filters frequency response
%I wish fdesign.bandpass to be reasonably close to the manual filter. Its clearly not!
x = 1: N;
freq = 0:(2*pi)/length(x):pi;
xdft = fft(output);
ydft = fft(maDiff);
figure;
plot(freq,abs(xdft(1:length(x)/2+1)),'m','linewidth',2);
hold on;
plot(freq,abs(ydft(1:length(x)/2+1)),'g','linewidth',2);
xlim([0 0.15]);
title(' Plot the filters frequency response');
legend('MA in cascade','fdesign.bandpass');

%%Plot the outputs in the time domain
%I wish fdesign.bandpass to be reasonably close to the manual filter. Its clearly not!
figure;
plot(data, 'r');
hold all;
plot(maDiff, 'm');  
plot(output, 'g');  
title(' Plot the filters time response');
legend('Raw Data','MA in cascade','fdesign.bandpass (Fc1)');

end