Solving problem of the Virtual bass enhancement (VBE) of audio

Hello recently I am trying to implement VBE for the enhancement of audio signal, I used the implementation of Non linear device(NLD) and Phase Vocoder(PV) . However the output seemS to be little bit distorted. Any one knows one of the best way to implement an NLD method with non distortion, or if he knows the implementation of maxxbass, Rbass or bass boost ? Please I need your help, thank you. Below is the matlab file for my NLD method. However the output sound tends to have some distortion like timbre has little bit changed. Anyone knows the best ways to implement a NLD?

%My VBE-NLD

clear all;

%close all

clc;

%%Load audio sample

%Fs=44100;

%audio=audioread('hiphop.wav');

%audio=audioread('sample.mp3');

%audio=audioread('audios/jazz.wav');

%audio=audioread('audios/classical.wav');

%audio=audioread('audios/hiphop.wav');

%audio=audioread('audios/rock.wav');

%[audio, Fs]=audioread('audios/drum.wav');

%[audio, Fs]=audioread('audios/8_micarray.wav');

[audio, Fs]=audioread('audios/song1.mp3');

%audio=audioread('audios/rock.wav');

%audio = sum(audio, 2) / size(audio, 2); %stereo to mono

%%%parameters

Fs=Fs;

x=audio;

Fcc=150; % cut-off frequency (Hz)

GdB=3; % NLD gain for transient calibration

%HPF

Fc=150; %cut-off frequency

N11=5;

[b,a]= High_pass(N11,Fc,Fs);

yt_high=Myiir_filter(b, a, x);

%LPF

Fc=150;

N11=5;

[b,a]= Low_pass(N11,Fc,Fs);

yt_low=Myiir_filter(b, a, x);

%y2=yt_high+yt_low;

%Apply gain

%G=10^(GdB/20);

%y3=G*yt_low;

%y4=yt_high+y3;

%y5=HPF_final(y4,Fc, Fs);

%NLD

H_hwr=[0.0278 0.5 1.8042 0 -6.0133 0 12.2658 0 -11.8891 0 4.3149];

H_fexp1=[0 1.4847 0 -1.4449 0 2.4713 0 -2.4234 0 0.9158 0];

h=H_hwr;

for n=1:11

if n==1

yt_low_proc=h(n);

else

yt_low_proc=yt_low_proc+h(n)*(yt_low.^(n-1));

end

end

h=H_fexp1;

for n=1:11

if n==1

yt_low_proc_2=h(n);

else

yt_low_proc_2=yt_low_proc_2+h(n)*(yt_low_proc.^(n-1));

end

end

%yt_low_proc_1=atsr(yt_low);

%yt_low_proc_2=abs(yt_low_proc_1);

%yt_low_proc_2=asqrt(yt_low);

%yt_low_proc_2=abs(yt_low_proc_1);

%yt_low_proc_2=my_full_wav_rect(yt_low);

%yt_low_proc_2=my_half_wav_rect(yt_low);

%my_softclipping(x, k)

%BPF

NBPF1 = 200; % Order

NBPF2 = 100; % Order

Fstop = 2*Fcc/4; % Stopband Frequency

Fpass = Fcc+ Fcc/4; % Passband Frequency

Fpass2=2*Fcc;

Fstop2=5*Fcc;

Wstop = 1; % Stopband Weight

Wstop2 = 0.001; % Stopband Weight

Wpass = 1; % Passband Weight

Wpass2 = 1; % Passband Weight

dens = 20; % Density Factor

b = firpm(NBPF1, [0 Fstop Fpass Fsreduced/2]/(Fsreduced/2), [0 0 1 1], [Wstop Wpass], ...

{dens});

Hd = dfilt.dffir(b);

yt_low_proc_bpf=filter(Hd,yt_low_proc_2);

b = firpm(NBPF2, [0 Fpass2 Fstop2 Fsreduced/2]/(Fsreduced/2), [1 1 0 0], [Wpass Wstop], ...

{dens});

Hd = dfilt.dffir(b);

yt_low_proc_bpf=filter(Hd,yt_low_proc_bpf);

%delay adjustment

delay_low=round((NBPF1+NBPF2)/2);

yt_high=[zeros(delay_low,1); transpose(yt_high)];

if length(yt_high)>length(yt_low_proc_bpf)

yt_high=yt_high(1:length(yt_low_proc_bpf));

else

yt_low_proc_bpf=yt_low_proc_bpf(1:length(yt_high));

end

%Apply gain

G=10^(GdB/20);

yt_low_proc_gain=G*yt_low_proc_bpf;

%Reconstruction

yt_proc=yt_high+transpose(yt_low_proc_gain);

%yt_proc=HPF_final(yt_proc, Fcc,Fs);

xfilt=MyHPF_final(x,Fcc,Fs);

yfilt=MyHPF_final(yt_proc,Fcc,Fs);

y_filt_low=MyLPF_final(x,Fcc,Fs);