Prev: A beginers question
Next: How to reduce the local illumination variations.
From: Wang on 13 Aug 2010 17:50
I've done with a implementation with a IIR Direct Form II Transposed filter, both in M-File code and Fortran 77 code. The algorithm is supposed to be completely the same as shown in the Matlab "Signal Processing Tool User's Guide". which is as follows

The M-File is as follows:
--------------------------------------------------------------
% Function: TDFII(b, a, x).
% Direct realizations of rational transfer functions.
% Input parameters:
% typ: 1 for direct realization, 2 for transposed
% b, a: numerator and denominator polynomials
% x: input sequence.
% Output:
% y: output sequence.
p = length(a)-1; q = length(b)-1; pq = max(p,q);
a = a(2:p+1); u = zeros(1,pq); % u: the internal state
for i = 1:length(x),
y(i) = b(1)*x(i)+u(1);
u = [u(2:pq),0];
u(1:q) = u(1:q) + b(2:q+1)*x(i);
u(1:p) = u(1:p) - a*y(i);
end

The Fortran 77 code is as follows:
--------------------------------------------------------
SUBROUTINE Filter( x, y, SL )
INTEGER SL ! Sample length
REAL x(SL), y(SL) ! x and y vector is input and output respectively

! Initialize the input
PARAMETER IIR_NSEC = 23
DOUBLE PRECISION in(SL), out(SL)
DOUBLE PRECISION CF_A(IIR_NSEC), CF_B(IIR_NSEC), z(IIR_NSEC-1)
INTEGER m, i, j, k, kk
SAVE z

!------------------------------------------------------
! Initialization of filter coefficients
!------------------------------------------------------
! Initializing the filter coeffcients, which should be constant while filtering
SAVE CF_A, CF_B
DATA CF_A
/ Here is the coefficeint generated by FDATools /

DATA CF_B
/ Here is the coefficeint generated by FDATools /

!---------------------------------------------------------
! Initialize the input and output vector for filtering
!---------------------------------------------------------
DO 5, i = 1, SL
in(i) = DBLE( x(i) )
out(i) = DBLE( 0 )
5 CONTINUE

!------------------------------------------------------
! filtering the signal vector
!------------------------------------------------------

! Reset past status of delay output
DO 10, j = 1, IIR_NSEC-1
z(j) = 0.0
10 CONTINUE

!---------------------------------
! Start filtering the k-th sample
! Omit CF_A(1) since it has been normalized to 1
DO 30, m = 1, SL
out(m) = CF_B(1) * in(m) + z(1)

! Update the past status of delay output
DO 20, kk=1, IIR_NSEC-2
z(kk) = z(kk+1)
20 CONTINUE
z(IIR_NSEC-1) = 0.0

! Calculate the delay outputs z(k)
DO 15, k = 1, IIR_NSEC-1
z(k) = z(k) + CF_B(k+1)*in(m) - CF_A(k+1)*out(m)
15 CONTINUE

30 CONTINUE

DO 50, im = 1, SL
y(SL) = REAL( out(SL) )
50 CONTINUE

END

I use the IDENTICAL coefficient both in M code and Fortran code, but the M code works fine in Matlab and the filtered signal becomes divergence after only about 19th loop.

I’ve tried the Fortran complier of Intel Visual Fortran 10 and Compaq Digital Fortran 6, almost the same error occurred during runtime.

Is it because the precision of data? I did found some difference between the results come from the same expression in M-code and Fortran code.

Please help me about the solution !!! This problem becomes a hurdle to my project! So somewhat rush!
 | 
Pages: 1
Prev: A beginers question
Next: How to reduce the local illumination variations.