butter_iir.m

DSP中巴特沃思滤波器的设计使用Verilog编写.

FilterOrder = 4;
data_width = 12;    % Data input bit width
max_data = 2^(data_width-1)-1;
min_data = -max_data;
biquad_number = floor((FilterOrder+1)/2);

% Filter Coefficients
%FeedForward = [0.2958    0.5876    0.2958];
%FeedBack = [1.0000    0.2138    0.4518];
[FeedForward,FeedBack] = butter(FilterOrder,300/500);
figure(1);
freqz(FeedForward,FeedBack);
title('Frequency Response (full precision coefficients)');
figure(2);
zplane(FeedForward,FeedBack);
title('Zero-pole plot using full precision coefficients');
[z,p,k] = tf2zp(FeedForward,FeedBack);

[SecOrd,G] = zp2sos(z,p,k);
G=0.16699


  
%Impulse response

   Inp = zeros(300,1);
   Inp(10) = max_data;
   %Impulse Response Plotting
   figure(3);
   biquad_output=zeros(biquad_number,min(length(Inp),300));
   biquad_output(1,:) = butter_biquad(SecOrd(1,:), Inp * 2^4);
   for i=1: biquad_number - 1;
       biquad_output(i+1,:) = butter_biquad(SecOrd(i+1,:), biquad_output(i,:));
   end
   output = biquad_output(biquad_number,:) / 2^4 * G;
   i_output = filter(FeedForward, FeedBack, Inp);
   impulse_output = output;
   sz = size(output);
   max_x = max(sz);
   xline = linspace (0, max_x - 1, max_x);
   plot (xline, output,'r');
   title ('Time Display of Impulse Response');
   grid on;
   zoom on;
   outfile1 = fopen('imp_in.txt','w');
	fprintf(outfile1, '%d\n', Inp);
   fclose(outfile1);
   outfile1 = fopen('m_Inpulse_out.txt','w');
	fprintf(outfile1, '%d\n', output);
   fclose(outfile1);
   
   % Frequency Domain Plotting 
    to_plot = output;
    figure(4);
    freqdat = fft(to_plot);
    absdat = abs(freqdat);
    maxdat = max (absdat);
    logdat = 20*log10(absdat);
    sz = size(to_plot);
    numpts = max(sz);
    freq_res = 1/numpts;
    xline = linspace (0, ( (1/2)- freq_res ), round(numpts/2) );
   	plot (xline, logdat(1: round(numpts/2) ), 'b');
    title ('Frequency Display of Impulse Response');
    grid on;
    zoom on;
    xlabel ('Frequency');
    ylabel ('Magnitude - dB');

    
%Step Response
 
    step= zeros(300,1);
 	for	i=10:300
        step(i) = max_data;
	end;
	%step Response Plotting
    figure(5);
    biquad_output=zeros(biquad_number,min(length(step),300));
    biquad_output(1,:) = butter_biquad(SecOrd(1,:), step * 2^4);
   for i=1: biquad_number - 1;
       biquad_output(i+1,:) = butter_biquad(SecOrd(i+1,:), biquad_output(i,:));
   end
   output = biquad_output(biquad_number,:) / 2^4 * G;
    s_output = filter(FeedForward, FeedBack, step);
    step_output = output;
    sz = size(output);
    max_x = max(sz);
    xline = linspace (0, max_x - 1, max_x);
  	plot (xline, output,'r');
    title ('Time Display of Step Response');
    grid on;
    zoom on;
    outfile2 = fopen('step_in.txt','w');
	fprintf(outfile2, '%d\n', step);
    fclose(outfile2);
      outfile2 = fopen('m_Step_out.txt','w');
	fprintf(outfile2, '%d\n', output);
    fclose(outfile2);
  
   % Frequency Domain Plotting 
    to_plot = output;
    figure(6);
    freqdat = fft(to_plot);
    absdat = abs(freqdat);
    maxdat = max (absdat);
    logdat = 20*log10(absdat);
    sz = size(to_plot);
    numpts = max(sz);
    freq_res = 1/numpts;
    xline = linspace (0, ( (1/2)- freq_res ), round(numpts/2) );
  	plot (xline, logdat(1: round(numpts/2) ), 'b');
    title ('Frequency Display of Step Response');
    grid on;
    zoom on;
    xlabel ('Frequency');
    ylabel ('Magnitude - dB');
    
%Ramdom Input

    random = round((rand(1, 100)-0.5)*(2^data_width-2)-1); 
    int_random = round(random);
    zero_init = zeros(10,1);
    zero_init = zero_init';
    int_random = [zero_init int_random];
    % Saturating Input Value
    a_in = find (int_random > max_data);
    b_in = find (int_random < min_data);
    if (~isempty(a_in)|~isempty(b_in))
         lenax = length (a_in);
         lenbx = length (b_in);
        for i = 1:lenax
             int_random(a_in(i)) = max_data;
        end
        for i = 1:lenbx
             int_random(b_in(i)) = min_data;
         end
    end
   %Random Response Plotting
   % This is the stimulus used in random.vwf
    int_random = [0 0 0 0 0 0 0 0 0 0 588 -1735 48 -461 920 1837 -930 -545 240 -1410 -420 518 193 -682 553 0 0 0 0 0 0];
    figure(7);
    biquad_output=zeros(biquad_number,min(length(int_random),300));
    biquad_output(1,:) = butter_biquad(SecOrd(1,:), int_random * 2^4);
   for i=1: biquad_number - 1;
       biquad_output(i+1,:) = butter_biquad(SecOrd(i+1,:), biquad_output(i,:));
   end
   output = biquad_output(biquad_number,:) / 2^4 * G;
   r_output = filter(FeedForward, FeedBack, int_random);
    random_output = output;
    sz = size(output);
    max_x = max(sz);
    xline = linspace (0, max_x - 1, max_x);
  	plot (xline, output,'r');
    title ('Time Display of random Response');
    grid on;
    zoom on;
    outfile3 = fopen('rand_in.txt','w');
	fprintf(outfile3, '%d\n', int_random);
	fclose(outfile3);
    outfile3 = fopen('m_Random_out.txt','w');
	fprintf(outfile3, '%d\n', output);
	fclose(outfile3);

  
    % Frequency Domain Plotting 
    to_plot = output;
    figure(8);
    freqdat = fft(to_plot);
    absdat = abs(freqdat);
    maxdat = max(absdat);
    logdat = 20*log10(absdat);
    sz = size(to_plot);
    numpts = max(sz);
    freq_res = 1/numpts;
    xline = linspace (0, ( (1/2)- freq_res ), round(numpts/2) );
   	plot (xline, logdat(1: round(numpts/2) ), 'b');
    title ('Frequency Display of Random Response');
    grid on;
    zoom on;
    xlabel ('Frequency');
    ylabel ('Magnitude - dB');