FIR design - convolution method

I came across the following problem in FIR design using window and convolution techniques .

Project a FIR filter using a Hamming window to achieve the following specifications:

Passband 0,3-3,4 kHz
Stopband 0-0,2 e 4-8 kHz
Stopband atenuation > 25 dB
Sampling frequency 32 kHz

Since the transition band is different , in order to solve the problem it is necessary to build a low pass filter, then a highpass filter and apply convolution to their coeficient, with this procedure it is possible to find a bandpass filter. Performing these steps, I found a correct magnitude plot but a weird phase plot.

Below I present the matlab code. I would appreciate if someone running it could realize what is wrong .

In spite of the phase found is linear as required by FIR filters, it is inverse from what it should be, besides , it is all over the band up to Fs/2, usually the phase is linear just in the passband, that’s why I think there is something wrong.

It is an academic problem I have to build a code and I am not allowed to use fir1 comand, neither remez , fir2 or any other , I Know that I can solve this problem using fir1 command, but that is not the point, the solution must be given through convolutioning the two filter coefficients, so , I must find out what is wrong with the matlab code wich follows:

%general data
Ap= 3;
As= 25;
Fs=32000;

% taking the smaller ripple
dp=(10^(Ap/20))-1; ds=10^(-As/20);
d=min(dp,ds);
As= -20*log10(d);

%Lowpass:

fp =3400;
fs= 4000;


% normalizing frequencies

fpn=fp/Fs;
fsn=fs/Fs;
wcn=2*pi*fpn;

% Filter order to Hamming window
N=ceil(3.3/(fsn-fpn));

% Hamming window
wd=hamming(N-1);

% Ideal filter
na=1N-1)/2;
hd=2*fpn*sin(na.*wcn)./(na.*wcn);
hdn=[fliplr(hd) 2*fpn hd];

%windowing
hn=hdn.*wd';


% Highpass

fp1=300;
fs1=200;

% normalizing frequencies

fpn1=fp1/Fs;
fsn1=fs1/Fs;
wcn1=2*pi*fpn1;

% Filter order to Hamming window
N1=ceil(3.3/(fpn1-fsn1));

% Hamming window
wd1=hamming(N1);

% Ideal filter
na1=1N1-1)/2;
hd1=-2*fpn1*sin(na1.*wcn1)./(na1.*wcn1);
hdn1=[fliplr(hd1) 1-2*fpn1 hd1];

%windowing
hn1=hdn1.*wd1';

% bandpass filter coeficients

H=conv(hn1,hn);

%Plotting frequency response

freqz(H,1,512,Fs);