Parks-McClellan method for best root-Nyquist filter search

rnyq_search.m

%Clear all previous variables etc.
clear; clc;

%The excess bandwidth for our reference square root raised cosine filter is
%set to alpha (a), number of taps is _span*sps+1_
%We set the gain of this filter to 1.0 at the passband
%Next, we calculate raised cosine filter based on our reference taps
%and compute the RMS ISI value
global a;
global span;
global sps;
a=0.5;
span=10;
sps=10;

%We divide by sqrt(sps) to obtain a unity gain filter
ref_rrc=rcosdesign(a, span, sps)/sqrt(sps);
ref_rc=conv(ref_rrc, ref_rrc);
ref_rc_isi_rms=rms([ref_rc(1:sps:span*sps),ref_rc(span*sps+sps+1:sps:2*span*sps+1)]);

%Place for optimal values
opt_vals=zeros(1,3);
init=[1.0, 1.0, 1.0];

%The real action takes place here
[opt_vals, value] = fminsearch(@rms_isi, init);

%We use Parks-McClellan method to obtain the filter coefficients
f = [0, opt_vals(2)*(1-a)*(1/sps), (1/sps), (1/sps), opt_vals(3)*(1+a)*(1/sps), 1];
am = [1, 1, 1/sqrt(2), 1/sqrt(2), 0, 0];
wv=[opt_vals(1),1,1];
b = firpm(span*sps,f,am,wv);

%Calculate Nyquist filter from root-Nyquist one
bb=conv(b,b);

%Freq response of both proposed and reference filters
fvtool(ref_rrc, 1, b, 1);

%Plot ISI for both proposed and reference filter
stem(bb(1:sps:2*span*sps+1), 'red o');
hold on;
stem(ref_rc(1:sps:2*span*sps+1), 'blue +');
ylim([-0.2e-3, +0.2e-3]);
xlim([0,2*span+1+1]);

%display RMS values for the proposed, reference and mismatched filters
mm=conv(b,ref_rrc);
fprintf('RMS ISI values:\n');
fprintf('Ref RC:\t\t\t%1.6f\n', ref_rc_isi_rms);
fprintf('Proposed RC:\t%1.6f\n',rms([bb(1:sps:span*sps), bb(span*sps+sps+1:sps:2*span*sps+1)]));
fprintf('Mismatched:\t\t%1.6f\n',rms([mm(1:sps:span*sps), mm(span*sps+sps+1:sps:2*span*sps+1)]));
fprintf('Gain over ref. RC: %1.1f\n\n', ref_rc_isi_rms/rms([bb(1:sps:span*sps), bb(span*sps+sps+1:sps:2*span*sps+1)]));
[h,~]=freqz(ref_rc,1000);
fprintf('Max passband ripple of ref. RC:\t%1.4f dB\n', max(20*log10(abs(h(1:100)*1000))));
[h,~]=freqz(bb,1000);
fprintf('Max passband ripple:\t\t\t%1.4f dB\n', max(20*log10(abs(h(1:100)*1000))));

%This function returns the root mean square value of the intersymbol
%interference values for the full width of the filter
%x(1) - first weight passed to Parks-McClellan method (for the passband)
%x(2) - beta_1 value, a parameter for the left transition band edge
%x(3) - beta_2 value, a parameter for the right transition band edge
function f = rms_isi(x)
    global a;
    global span;
    global sps;
    f = [0, x(2)*(1-a)*(1/sps), (1/sps), (1/sps), x(3)*(1+a)*(1/sps), 1];
    am = [1, 1, 1/sqrt(2), 1/sqrt(2), 0, 0];
    wv=[x(1),1,1];
    
    %Calculate filter coefficients
    b = firpm(span*sps,f,am,wv);
    
    %Find ISI RMS
    bb=conv(b,b);
    f=rms([bb(1:sps:span*sps),bb(span*sps+sps+1:sps:2*span*sps+1)]);
end