Code covered by the BSD License  

Highlights from
Numerically Design Power Complementary IIR(FIR) Filter

Numerically Design Power Complementary IIR(FIR) Filter



02 May 2012 (Updated )

One iir filter is given, design power (or magnitude) complementary filter numerically

% This script finds the amplitude complementary IIR filter.
% 1st IIR filter is designed using matlab function 2nd IIR filter is found
% by using error minimization using fminsearch matlab function.

format long g;
Fs                    = 44100;       % sampling frequency of the underlying signal in Hz
freqPoints            = 1024;

% Design maximally flat Lowpass filter using matlab command
N                     = 4;         % order of the filter
Fc                    = 10000;        % cut off frequency in Hz
[b,a,b1,b2,sos_var,g] = maxflat(N,N,Fc/(Fs/2));
%[b,a,b1,b2,sos_var,g]= maxflat(N, Fc/(Fs/2));
Hd                    = dfilt.df2sos(sos_var, g);
[b, a]                = tf(Hd);      % b:numerator coefficient, a:denominator coefficient

% Design Magnitude Complementary filter using fminsearch matlab command

% put num and den in a single input vector for Lowpass filter
coeffLow              = [b,a];       

% intitlize 2nd (Highpass)filter with reversed and randomized 1st filter's coefficients.
coeffHigh             = [-b(end:-1:1)+0.01*rand(1,length(b)),-a(end:-1:1)+0.01*rand(1,length(a))];

% find 2nd filter's coefficients such that  filer coeffHigh and coeffLow minimises a cost function.
[coeffHigh, fval, exitflag] = fminsearch(@(coeffHigh)costFunctionIIR(coeffLow,coeffHigh,freqPoints), coeffHigh);

% run some iterations...
for i=1:5,
    %[y,fval, exitflag] = fminsearch(@(y)costFunctionIIR(b,y,freqPoints), x);
    [coeffHigh, fval, exitflag] = fminsearch(@(coeffHigh)costFunctionIIR(coeffLow,coeffHigh,freqPoints), t);

% Testing the results.
m0=(abs(freqz(coeffLow(1:end/2),coeffLow(end/2+1:end),freqPoints))); % 1st filter magnitude response
m1=(abs(freqz(coeffHigh(1:end/2),coeffHigh(end/2+1:end),freqPoints))); % 2nd filter magnitude response
plot(w,m0);hold on;  % plot 1st filter
plot(w,m1,'r'); % plot 2nd filter
plot(w,m0.*m0 + m1.*m1,'g'); % plot summed responce
hold off
beep;beep % signal after comlpetion.

Contact us