pfir_design_example.m

clear all;
FilterOrder = 128-1; % Make odd to force linear phase
Bands = 2;
SampleRate = 1e6;
N = 2^20;

% Some filter we want to compensate for
b = 2^(-14)*[-53 0 313 0 -1155 0 4989 8192 4989 0 -1155 0 313 0 -53];
a = 1;

% Create filter response
[response,frequencies] = freqz(b,a,N,SampleRate);
response = abs(response);

frequenciesNorm = frequencies./(SampleRate/2);

% Create inverse target up to a certain frequency
passPercent = 0.6;
stopIndx = floor(N*passPercent);
A1 = 1./response(1:stopIndx);
F1 = frequenciesNorm(1:stopIndx);
A2 = zeros(N-stopIndx,1);
F2 = frequenciesNorm(stopIndx+1:N);

% Set weight for what we favor more
W1 = 0.1      .*ones(stopIndx,1).';
W2 = 0.00001  .*ones(1,length(A2));

% Design filter
d = fdesign.arbmag('N,B,F,A',FilterOrder,Bands,F1,A1,F2,A2);
Hd = design(d,'equiripple','B1Weights',W1,'B2Weights',W2,'SystemObject',true);

% Calculate new filter's response
[responseComp, frequencies] = freqz(Hd,N,SampleRate);

%% Combine responses and view both
combined = 20*log10(responseComp.*response);
response = 20*log10(response);
plot(frequencies, combined, frequencies, response);
ylim([ -50 5 ]);
legend('Combined', 'Original');
xlabel('Frequency (Hz)'); ylabel('Magnitude (dB)');
grid on;