Software Implementation of LTI System from its Constant-Coefficient Difference Equation

A Linear time-invariant system is a system that obeys the following fundamental properties:

  1. Linearity: The input and output relationship of the system is based on linear differenctial eqautions. If the system maps inputs x1(t) and x2(t) to outputs y1(t) and y2(t) respectively, then it will map x3(t) = x1(t) + x2(t) to the output y3(t) where y3(t) = y1(t) + y2(t). (superposition principle)
  2. Time invariance: Time invariance means that the output of the system does not depend on the particular time the input is applied.

The fundamental result of the LTI system is that it can be characterized entirely by the systems impulse response.

The output of the system can be obtained by convolving the input with the impulse response or vice versa.

This article aims at developing a simple software library that can realise any LTI system from its Constant-Coefficient Difference Equation defined as:

3

The above equation can be represented using the following diagram (Direct 1 form):

3

Non-Recursive System,

3

Recursive System,

3

If the recursive and non-recursive part of the above general Constant-Coefficient Difference Equation is interchanged the overall system response remains the same as the convolution of two impulse response is commutative.

Thus the direct 1 form can be changed into direct 2 form:

3

Thus, recursive System,

3

Non-Recursive System,

3

The above form reduces the number of delays (memory) required to implement the system. (No. of delays = max{M, N} delays.)

The following code snippet implements the above form of LTI System that can be configured depending upon the difference equation of the system.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
 
#include <stdlib.h>

#include <math.h>
#include <stdio.h>

typedef struct _lti_sys {

    /* data */
    int M, N, delay_mem_size;
    float *coeff_buffer_a, *coeff_buffer_b;
    float *__delay_mem_buffer;
    int __delay_buffer_head_index;

} LTI_Sys;

int LTI_Sys_Init(LTI_Sys *lti_sys_handle, int M, int N);
int LTI_Sys_Set_Initial_State(LTI_Sys *lti_sys_handle, float *init_state, int size);
int LTI_Sys_Set_Coeff(LTI_Sys *lti_sys_handle, float *coeff_a, float *coeff_b, int size_a, int size_b);
int LTI_Sys_Process(LTI_Sys *lti_sys_handle, float ip, float *op);
int LTI_Sys_Reset(LTI_Sys *lti_sys_handle);
int LTI_Sys_Free(LTI_Sys *lti_sys_handle);

static int __get_delay_mem(LTI_Sys *lti_sys_handle, int index, float *op_val);
static int __update_delay_mem(LTI_Sys *lti_sys_handle, float data_value);

int LTI_Sys_Init(LTI_Sys *lti_sys_handle, int M, int N) {

    if (M == 0 && N == 0)
        return -1;

    lti_sys_handle->delay_mem_size = M >= N ? M : N;
    lti_sys_handle->M = M + 1;
    lti_sys_handle->N = N;
    lti_sys_handle->__delay_buffer_head_index = 0;
    lti_sys_handle->__delay_mem_buffer = NULL;
    lti_sys_handle->coeff_buffer_a = NULL;
    lti_sys_handle->coeff_buffer_b = NULL;

    if((lti_sys_handle->__delay_mem_buffer = malloc(sizeof(float) * lti_sys_handle->delay_mem_size)) == NULL) 
        goto free_n_exit_error;

    for(int i = 0; i < lti_sys_handle->delay_mem_size; ++i) {
        lti_sys_handle->__delay_mem_buffer[i] = 0.;
    }

    if(lti_sys_handle->N > 0) {
        if((lti_sys_handle->coeff_buffer_a = malloc(sizeof(float) * lti_sys_handle->N)) == NULL)
            goto free_n_exit_error;
        
        for(int i = 0; i < lti_sys_handle->N; ++i) {
            lti_sys_handle->coeff_buffer_a[i] = 0.;
        }
    }

    if(lti_sys_handle->M > 0) {
        if((lti_sys_handle->coeff_buffer_b = malloc(sizeof(float) * lti_sys_handle->M)) == NULL)
            goto free_n_exit_error;

        for(int i = 0; i < lti_sys_handle->M; ++i) {
            lti_sys_handle->coeff_buffer_b[i] = 0.;
        }
    }

    return 0;

    free_n_exit_error:
        free(lti_sys_handle->__delay_mem_buffer);
        free(lti_sys_handle->coeff_buffer_a);
        free(lti_sys_handle->coeff_buffer_b);
        return -1;

}

int LTI_Sys_Set_Initial_State(LTI_Sys *lti_sys_handle, float *init_state, int size) {

    if(size > lti_sys_handle->delay_mem_size)
        return -1;

    for(int i = 0; i < size; ++i) {
        lti_sys_handle->__delay_mem_buffer[lti_sys_handle->delay_mem_size - i - 1] = init_state[i];
    }

    return 0;

}

int LTI_Sys_Set_Coeff(LTI_Sys *lti_sys_handle, float *coeff_a, float *coeff_b, int size_a, int size_b) {

    if(size_a > lti_sys_handle->N || size_b > lti_sys_handle->M)
        return -1;

    for(int i = 0; i < size_a; ++i)
        lti_sys_handle->coeff_buffer_a[i] = coeff_a[i];
    for(int i = 0; i < size_b; ++i)
        lti_sys_handle->coeff_buffer_b[i] = coeff_b[i];

    return 0;

}

int LTI_Sys_Process(LTI_Sys *lti_sys_handle, float ip, float *op) {

    float op_val = ip, tmp_coeff = 0., tmp_op_val = 0.;

    for(int i = 0; i < lti_sys_handle->N; ++i) {
        if(!__get_delay_mem(lti_sys_handle, i, &tmp_coeff))
            op_val += lti_sys_handle->coeff_buffer_a[i] * tmp_coeff;
        else
            return -1;
    }

    tmp_op_val = op_val;

    for(int i = 0; i < lti_sys_handle->M; ++i) {
        if(i == 0)
            op_val = lti_sys_handle->coeff_buffer_b[i] * op_val;
        else if(!__get_delay_mem(lti_sys_handle, i - 1, &tmp_coeff))
            op_val += lti_sys_handle->coeff_buffer_b[i] * tmp_coeff;
        else
            return -1;
    }

    if(__update_delay_mem(lti_sys_handle, tmp_op_val))
        return -1;

    *op = op_val;

    return 0;

}

int LTI_Sys_Reset(LTI_Sys *lti_sys_handle) {

    lti_sys_handle->__delay_buffer_head_index = 0;

    for(int i = 0; i < lti_sys_handle->delay_mem_size; ++i) {
        lti_sys_handle->__delay_mem_buffer[i] = 0.;
    }

    for(int i = 0; i < lti_sys_handle->N; ++i) {
        lti_sys_handle->coeff_buffer_a[i] = 0.;
    }

    for(int i = 0; i < lti_sys_handle->M; ++i) {
        lti_sys_handle->coeff_buffer_b[i] = 0.;
    }

    return 0;

}

int LTI_Sys_Free(LTI_Sys *lti_sys_handle) {

    LTI_Sys_Reset(lti_sys_handle);

    free(lti_sys_handle->__delay_mem_buffer);
    free(lti_sys_handle->coeff_buffer_a);
    free(lti_sys_handle->coeff_buffer_b);
    return 0;

}

static int __get_delay_mem(LTI_Sys *lti_sys_handle, int index, float *op_val) {

    *op_val = lti_sys_handle->__delay_mem_buffer[(index + lti_sys_handle->__delay_buffer_head_index) % lti_sys_handle->delay_mem_size];
    return 0;

}

static int __update_delay_mem(LTI_Sys *lti_sys_handle, float data_value) {

    if(lti_sys_handle->__delay_buffer_head_index == 0)
        lti_sys_handle->__delay_buffer_head_index = lti_sys_handle->delay_mem_size - 1;
    else
        lti_sys_handle->__delay_buffer_head_index -= 1;

    lti_sys_handle->__delay_mem_buffer[lti_sys_handle->__delay_buffer_head_index] = data_value;
    return 0;

}

Example

Example Implementation of moving average filter: lti.c, lti.h

References