dsp.h

Go to the documentation of this file.

/*   libDSP - a digital signal processing library
 *   Copyright © 2017-2022  Ilia Platone
 *
 *   This program is free software; you can redistribute it and/or
 *   modify it under the terms of the GNU Lesser General Public
 *   License as published by the Free Software Foundation; either
 *   version 3 of the License, or (at your option) any later version.
 *
 *   This program is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *   Lesser General Public License for more details.
 *
 *   You should have received a copy of the GNU Lesser General Public License
 *   along with this program; if not, write to the Free Software Foundation,
 *   Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 */
 
#ifndef _DSP_H
#define _DSP_H
 
#ifdef  __cplusplus
extern "C" {
#endif
 
#ifndef DLL_EXPORT
#ifdef _WIN32
#define DLL_EXPORT __declspec(dllexport)
#else
#define DLL_EXPORT extern
#endif
#endif
 
#if defined(__linux__) || defined(_WIN32) || defined(__CYGWIN__)
#include <endian.h>
#else
#define __bswap_16(a) __builtin_bswap16(a)
#define __bswap_32(a) __builtin_bswap32(a)
#define __bswap_64(a) __builtin_bswap64(a)
#endif
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include <sys/types.h>
#include <time.h>
#include <assert.h>
#include <pthread.h>
#include <time.h>
#ifndef M_PI
    #define M_PI 3.14159265358979323846
#endif
 
#define DSP_MAX_STARS 200
typedef double dsp_t;
typedef double complex_t[2];
#define dsp_t_max 255
#define dsp_t_min -dsp_t_max
#define DSP_NAME_SIZE 128
 
DLL_EXPORT unsigned long int dsp_max_threads(unsigned long value);
 
#ifndef DSP_DEBUG
#define DSP_DEBUG
DLL_EXPORT void dsp_set_debug_level(int value);
DLL_EXPORT int dsp_get_debug_level();
DLL_EXPORT void dsp_set_app_name(char* name);
DLL_EXPORT char* dsp_get_app_name();
DLL_EXPORT void dsp_set_stdout(FILE *f);
DLL_EXPORT void dsp_set_stderr(FILE *f);
 
DLL_EXPORT void dsp_print(int x, char* str);
 
#define DSP_DEBUG_INFO 0
#define DSP_DEBUG_ERROR 1
#define DSP_DEBUG_WARNING 2
#define DSP_DEBUG_DEBUG 3
#define pdbg(x, ...) ({ \
char str[500]; \
struct timespec ts; \
time_t t = time(NULL); \
struct tm tm = *localtime(&t); \
clock_gettime(CLOCK_REALTIME, &ts); \
sprintf(str, "[%04d-%02d-%02dT%02d:%02d:%02d.%03ld ", tm.tm_year + 1900, tm.tm_mon + 1, tm.tm_mday, tm.tm_hour, tm.tm_min, tm.tm_sec, ts.tv_nsec/1000000); \
switch(x) { \
    case DSP_DEBUG_ERROR: \
    sprintf(&str[strlen(str)], "ERRO]"); \
        break; \
    case DSP_DEBUG_WARNING: \
    sprintf(&str[strlen(str)], "WARN]"); \
        break; \
    case DSP_DEBUG_DEBUG: \
    sprintf(&str[strlen(str)], "DEBG]"); \
        break; \
    default: \
    sprintf(&str[strlen(str)], "INFO]"); \
        break; \
} \
if(dsp_get_app_name() != NULL) \
    sprintf(&str[strlen(str)], "[%s]", dsp_get_app_name()); \
sprintf(&str[strlen(str)], " "); \
sprintf(&str[strlen(str)], __VA_ARGS__); \
dsp_print(x, str); \
})
#define pinfo(...) pdbg(DSP_DEBUG_INFO, __VA_ARGS__)
#define perr(...) pdbg(DSP_DEBUG_ERROR, __VA_ARGS__)
#define pwarn(...) pdbg(DSP_DEBUG_WARNING, __VA_ARGS__)
#define pgarb(...) pdbg(DSP_DEBUG_DEBUG, __VA_ARGS__)
#define pfunc pgarb("%s\n", __func__)
#define start_gettime
#define end_gettime
#else
#define pinfo(...)
#define perr(...)
#define pwarn(...)
#define pgarb(...)
#define pfunc(...)
#define start_gettime(...)
#define end_gettime(...)
#endif
 
 
#ifndef Min
#define Min(a,b) \
   ({ __typeof (a) _a = (a); \
       __typeof (a) _b = (b); \
     _a < _b ? _a : _b; })
#endif
#ifndef Max
#define Max(a,b) \
   ({ __typeof (a) _a = (a); \
       __typeof (a) _b = (b); \
     _a > _b ? _a : _b; })
#endif
#ifndef Log
#define Log(a,b) \
( log(a) / log(b) )
#endif
#ifndef DSP_ALIGN_TRANSLATED
#define DSP_ALIGN_TRANSLATED 1
#endif
#ifndef DSP_ALIGN_SCALED
#define DSP_ALIGN_SCALED 2
#endif
#ifndef DSP_ALIGN_ROTATED
#define DSP_ALIGN_ROTATED 4
#endif
#ifndef DSP_ALIGN_NO_MATCH
#define DSP_ALIGN_NO_MATCH 8
#endif
typedef struct dsp_point_t
{
    double* location;
    int dims;
} dsp_point;
 
typedef struct dsp_offset_t
{
    double* offset;
    int dims;
} dsp_offset;
 
typedef struct dsp_star_t
{
    dsp_point center;
    double diameter;
    double peak;
    double flux;
    double theta;
    char name[DSP_NAME_SIZE];
} dsp_star;
 
typedef struct dsp_triangle_t
{
    double index;
    int dims;
    double *theta;
    double *sizes;
    double *ratios;
    dsp_star *stars;
    int stars_count;
} dsp_triangle;
 
typedef struct dsp_align_info_t
{
    double* offset;
    double* center;
    double* radians;
    double* factor;
    int dims;
    dsp_triangle triangles[2];
    int triangles_count;
    double score;
    double decimals;
    int err;
} dsp_align_info;
 
typedef union
{
    struct
    {
        double real;
        double imaginary;
    } *complex;
    complex_t *pairs;
    double *buf;
} dsp_complex;
 
typedef struct dsp_region_t
{
    int start;
    int len;
} dsp_region;
 
typedef union dsp_location_t
{
    struct
    {
        double x;
        double y;
        double z;
    } xyz;
    struct
    {
        double lon;
        double lat;
        double el;
    } geographic;
    double coordinates[3];
} dsp_location;
 
typedef void *(*dsp_func_t) (void *, ...);
 
typedef struct dsp_stream_t
{
    char name[DSP_NAME_SIZE];
    int is_copy;
    int len;
    int dims;
    int* sizes;
    dsp_t* buf;
    dsp_complex dft;
    void *arg;
    struct dsp_stream_t* parent;
    struct dsp_stream_t** children;
    int child_count;
    dsp_location* location;
    double* target;
    struct timespec starttimeutc;
    double wavelength;
    double focal_ratio;
    double diameter;
    double SNR;
    int red;
    double *pixel_sizes;
    double samplerate;
    pthread_t thread;
    dsp_func_t func;
    struct dsp_stream_t *magnitude;
    struct dsp_stream_t *phase;
    dsp_region *ROI;
    dsp_star *stars;
    int stars_count;
    dsp_triangle *triangles;
    int triangles_count;
    dsp_align_info align_info;
    int frame_number;
} dsp_stream, *dsp_stream_p;
 
DLL_EXPORT void dsp_fourier_dft(dsp_stream_p stream, int exp);
 
DLL_EXPORT void dsp_fourier_idft(dsp_stream_p stream);
 
DLL_EXPORT void dsp_fourier_2dsp(dsp_stream_p stream);
 
DLL_EXPORT void dsp_fourier_2complex_t(dsp_stream_p stream);
 
DLL_EXPORT void dsp_fourier_phase_mag_array_get_complex(double* mag, double* phi, complex_t *out, int len);
 
DLL_EXPORT double* dsp_fourier_complex_array_get_magnitude(dsp_complex in, int len);
 
DLL_EXPORT double* dsp_fourier_complex_array_get_phase(dsp_complex in, int len);
 
DLL_EXPORT void dsp_filter_squarelaw(dsp_stream_p stream);
 
DLL_EXPORT void dsp_filter_lowpass(dsp_stream_p stream, double frequency);
 
DLL_EXPORT void dsp_filter_highpass(dsp_stream_p stream, double frequency);
 
DLL_EXPORT void dsp_filter_bandpass(dsp_stream_p stream, double LowFrequency,
                                    double HighFrequency);
 
DLL_EXPORT void dsp_filter_bandreject(dsp_stream_p stream, double LowFrequency,
                                      double HighFrequency);
 
DLL_EXPORT void dsp_convolution_convolution(dsp_stream_p stream, dsp_stream_p matrix);
 
DLL_EXPORT void dsp_convolution_correlation(dsp_stream_p stream, dsp_stream_p matrix);
 
#ifndef dsp_stats_min
#define dsp_stats_min(buf, len)\
({\
    int i;\
    __typeof(buf[0]) min = (__typeof(buf[0]))buf[0];\
    for(i = 0; i < len; i++) {\
        min = Min(buf[i], min);\
    }\
    min;\
    })
#endif
 
#ifndef dsp_stats_max
#define dsp_stats_max(buf, len)\
({\
    int i;\
    __typeof(buf[0]) max = (__typeof(buf[0]))buf[0];\
    for(i = 0; i < len; i++) {\
        max = Max(buf[i], max);\
    }\
    max;\
    })
#endif
 
#ifndef dsp_stats_mid
#define dsp_stats_mid(buf, len)\
({\
    int i;\
    __typeof(buf[0]) min = dsp_stats_min(buf, len);\
    (__typeof(buf[0]))(min - dsp_stats_max(buf, len)) / 2.0 + min;\
})
#endif
 
#ifndef dsp_stats_maximum_index
#define dsp_stats_minimum_index(buf, len)\
({\
    int i;\
    __typeof(buf[0]) min = dsp_stats_min(buf, len);\
    for(i = 0; i < len; i++) {\
        if(buf[i] == min) break;\
    }\
    i;\
    })
#endif
 
#ifndef dsp_stats_maximum_index
#define dsp_stats_maximum_index(buf, len)\
({\
    int i;\
    __typeof(buf[0]) max = dsp_stats_max(buf, len);\
    for(i = 0; i < len; i++) {\
        if(buf[i] == max) break;\
    }\
    i;\
    })
#endif
 
#ifndef dsp_stats_stddev
#define dsp_stats_mean(buf, len)\
({\
    int __dsp__i;\
    double __dsp__mean = 0;\
    for(__dsp__i = 0; __dsp__i < len; __dsp__i++) {\
        __dsp__mean += buf[__dsp__i];\
    }\
    __dsp__mean /= len;\
    __dsp__mean;\
    })
#endif
 
#ifndef dsp_stats_stddev
#define dsp_stats_stddev(buf, len)\
({\
    double __dsp__mean = dsp_stats_mean(buf, len);\
    int __dsp__x;\
    double __dsp__stddev = 0;\
    for(__dsp__x = 0; __dsp__x < len; __dsp__x++) {\
        __dsp__stddev += fabs(buf[__dsp__x] - __dsp__mean);\
    }\
    __dsp__stddev /= len;\
    __dsp__stddev;\
    })
#endif
 
#ifndef dsp_stats_val_count
#define dsp_stats_val_count(buf, len, val) \
({\
    int x;\
    int count = 0;\
    for(x = 0; x < len; x++) {\
        if(buf[x] == val)\
            count ++;\
    }\
    count;\
    })
#endif
 
#ifndef dsp_stats_val_sum
#define dsp_stats_val_sum(buf, len) \
({\
    int x;\
    double sum = 0;\
    for(x = 0; x < len; x++) {\
        sum += buf[x];\
    }\
    sum;\
    })
#endif
 
#ifndef dsp_stats_range_count
#define dsp_stats_range_count(buf, len, lo, hi) \
({\
    int x;\
    int count = 0;\
    for(x = 0; x < len; x++) {\
        if(buf[x] < hi && buf[x] >= lo)\
            count ++;\
    }\
    count;\
    })
#endif
 
#ifndef dsp_stats_compare
#define dsp_stats_compare(in1, in2, len)\
({\
    __typeof(in1[0]) out = 0;\
    for(int i = 0; i < len; i++) {\
        out += in1[i] - (__typeof(in1[0]))in2[i];\
    }\
    out;\
    })
#endif
 
DLL_EXPORT double* dsp_stats_histogram(dsp_stream_p stream, int size);
 
DLL_EXPORT void dsp_buffer_shift(dsp_stream_p stream);
 
DLL_EXPORT void dsp_buffer_removemean(dsp_stream_p stream);
 
#ifndef dsp_buffer_stretch
#define dsp_buffer_stretch(buf, len, _mn, _mx)\
({\
    int k;\
    __typeof(buf[0]) __mn = dsp_stats_min(buf, len);\
    __typeof(buf[0]) __mx = dsp_stats_max(buf, len);\
    double oratio = (_mx - _mn);\
    double iratio = (__mx - __mn);\
    if(iratio == 0) iratio = 1;\
    for(k = 0; k < len; k++) {\
        buf[k] -= __mn;\
        buf[k] = (__typeof(buf[0]))((double)buf[k] * oratio / iratio);\
        buf[k] += _mn;\
    }\
})
#endif
 
#ifndef dsp_buffer_set
#define dsp_buffer_set(buf, len, _val)\
({\
    int k;\
    for(k = 0; k < len; k++) {\
        buf[k] = (__typeof(buf[0]))(_val);\
    }\
})
#endif
 
#ifndef dsp_buffer_normalize
#define dsp_buffer_normalize(buf, len, mn, mx)\
({\
    int k;\
    for(k = 0; k < len; k++) {\
        buf[k] = Max(mn, Min(mx, buf[k]));\
    }\
})
#endif
 
DLL_EXPORT void dsp_buffer_max(dsp_stream_p stream, dsp_t* in, int len);
 
DLL_EXPORT void dsp_buffer_min(dsp_stream_p stream, dsp_t* in, int len);
 
DLL_EXPORT void dsp_buffer_sub(dsp_stream_p stream, dsp_t* in, int len);
 
DLL_EXPORT void dsp_buffer_sum(dsp_stream_p stream, dsp_t* in, int len);
 
DLL_EXPORT void dsp_buffer_div(dsp_stream_p stream, dsp_t* in, int len);
 
DLL_EXPORT void dsp_buffer_mul(dsp_stream_p stream, dsp_t* in, int len);
 
DLL_EXPORT void dsp_buffer_pow(dsp_stream_p stream, dsp_t* in, int len);
 
DLL_EXPORT void dsp_buffer_log(dsp_stream_p stream, dsp_t* in, int len);
 
DLL_EXPORT void dsp_buffer_sub1(dsp_stream_p stream, dsp_t val);
 
DLL_EXPORT void dsp_buffer_1sub(dsp_stream_p stream, dsp_t val);
 
DLL_EXPORT void dsp_buffer_sum1(dsp_stream_p stream, dsp_t val);
 
DLL_EXPORT void dsp_buffer_div1(dsp_stream_p stream, double val);
 
DLL_EXPORT void dsp_buffer_1div(dsp_stream_p stream, double val);
 
DLL_EXPORT void dsp_buffer_mul1(dsp_stream_p stream, double val);
 
DLL_EXPORT void dsp_buffer_pow1(dsp_stream_p stream, double val);
 
DLL_EXPORT void dsp_buffer_log1(dsp_stream_p stream, double val);
 
DLL_EXPORT void dsp_buffer_median(dsp_stream_p stream, int size, int median);
 
DLL_EXPORT void dsp_buffer_sigma(dsp_stream_p stream, int size);
 
DLL_EXPORT void dsp_buffer_deviate(dsp_stream_p stream, dsp_t* deviation, dsp_t mindeviation, dsp_t maxdeviation);
 
#ifndef dsp_buffer_reverse
#define dsp_buffer_reverse(buf, len) \
    ({ \
        int i = (len - 1) / 2; \
        int j = i + 1; \
        __typeof(buf[0]) _x; \
        while(i >= 0) \
        { \
          _x = buf[j]; \
          buf[j] = buf[i]; \
          buf[i] = _x; \
          i--; \
          j++; \
        } \
    })
#endif
 
#ifndef dsp_buffer_swap
#define dsp_buffer_swap(in, len) \
    ({ \
        int k; \
        switch(sizeof(((__typeof (in[0])*)in)[0])) { \
        case 2: \
            for(k = 0; k < len; k++) \
                ((__typeof (in[0])*)in)[k] = __bswap_16(((__typeof (in[0])*)in)[k]); \
            break; \
        case 3: \
            for(k = 0; k < len; k++) \
            ((__typeof (in[0])*)in)[k] = __bswap_32(((__typeof (in[0])*)in)[k]); \
            break; \
        } \
    })
#endif
 
#ifndef dsp_buffer_copy
#define dsp_buffer_copy(in, out, len) \
    ({ \
        int k; \
        for(k = 0; k < len; k++) { \
        ((__typeof (out[0])*)out)[k] = (__typeof (out[0]))((__typeof (in[0])*)in)[k]; \
        } \
    })
#endif
 
#ifndef dsp_buffer_copy_stepping
#define dsp_buffer_copy_stepping(in, out, inlen, outlen, instep, outstep) \
    ({ \
    int k; \
    int t; \
        for(k = 0, t = 0; k < inlen && t < outlen; k+=instep, t+=outstep) { \
        ((__typeof (out[0])*)out)[t] = (__typeof (out[0]))((__typeof (in[0])*)in)[k]; \
        } \
    })
#endif
 
DLL_EXPORT void dsp_stream_alloc_buffer(dsp_stream_p stream, int len);
 
DLL_EXPORT void dsp_stream_set_buffer(dsp_stream_p stream, void *buffer, int len);
 
DLL_EXPORT dsp_t* dsp_stream_get_buffer(dsp_stream_p stream);
 
DLL_EXPORT void dsp_stream_free_buffer(dsp_stream_p stream);
 
DLL_EXPORT dsp_stream_p dsp_stream_new(void);
 
DLL_EXPORT void dsp_stream_free(dsp_stream_p stream);
 
DLL_EXPORT dsp_stream_p dsp_stream_copy(dsp_stream_p stream);
 
DLL_EXPORT void dsp_stream_add_child(dsp_stream_p stream, dsp_stream_p child);
 
DLL_EXPORT void dsp_stream_add_star(dsp_stream_p stream, dsp_star star);
 
DLL_EXPORT void dsp_stream_del_star(dsp_stream_p stream, int n);
 
DLL_EXPORT void dsp_stream_add_triangle(dsp_stream_p stream, dsp_triangle triangle);
 
DLL_EXPORT void dsp_stream_del_triangle(dsp_stream_p stream, int index);
 
DLL_EXPORT void dsp_stream_calc_triangles(dsp_stream_p stream);
 
DLL_EXPORT void dsp_stream_del_child(dsp_stream_p stream, int n);
 
DLL_EXPORT void dsp_stream_add_dim(dsp_stream_p stream, int len);
 
DLL_EXPORT void dsp_stream_set_dim(dsp_stream_p stream, int dim, int size);
 
DLL_EXPORT void dsp_stream_del_dim(dsp_stream_p stream, int n);
 
DLL_EXPORT int dsp_stream_set_position(dsp_stream_p stream, int *pos);
 
DLL_EXPORT int* dsp_stream_get_position(dsp_stream_p stream, int index);
 
DLL_EXPORT void *dsp_stream_exec(dsp_stream_p stream, void *args, ...);
 
DLL_EXPORT void dsp_stream_crop(dsp_stream_p stream);
 
DLL_EXPORT void dsp_stream_rotate(dsp_stream_p stream);
 
DLL_EXPORT void dsp_stream_sum(dsp_stream_p stream, dsp_stream_p in);
 
DLL_EXPORT void dsp_stream_multiply(dsp_stream_p stream, dsp_stream_p in);
 
DLL_EXPORT void dsp_stream_subtract(dsp_stream_p stream, dsp_stream_p in);
 
DLL_EXPORT void dsp_stream_translate(dsp_stream_p stream);
 
DLL_EXPORT void dsp_stream_scale(dsp_stream_p stream);
 
DLL_EXPORT void dsp_stream_align(dsp_stream_p stream);
 
DLL_EXPORT void dsp_signals_whitenoise(dsp_stream_p stream);
 
DLL_EXPORT void dsp_signals_sinewave(dsp_stream_p stream, double samplefreq, double freq);
 
DLL_EXPORT void dsp_signals_sawtoothwave(dsp_stream_p stream, double samplefreq, double freq);
 
DLL_EXPORT void dsp_signals_triwave(dsp_stream_p stream, double samplefreq, double freq);
 
DLL_EXPORT void dsp_modulation_frequency(dsp_stream_p stream, double samplefreq, double freq, double bandwidth);
 
DLL_EXPORT void dsp_modulation_amplitude(dsp_stream_p stream, double samplefreq, double freq);
 
void dsp_recons_align(dsp_stream_p stream, dsp_stream_p matrix);
 
DLL_EXPORT dsp_stream_p* dsp_file_read_fits(const char* filename, int *channels, int stretch);
 
DLL_EXPORT void dsp_file_write_fits(const char* filename, int bpp, dsp_stream_p stream);
 
DLL_EXPORT void dsp_file_write_fits_composite(const char* filename, int components, int bpp, dsp_stream_p* stream);
 
DLL_EXPORT dsp_stream_p* dsp_file_read_jpeg(const char* filename, int *channels, int stretch);
 
DLL_EXPORT void dsp_file_write_jpeg(const char* filename, int quality, dsp_stream_p stream);
 
DLL_EXPORT void dsp_file_write_jpeg_composite(const char* filename, int components, int quality, dsp_stream_p* stream);
 
DLL_EXPORT dsp_stream_p* dsp_file_read_png(const char* filename, int *channels, int stretch);
 
DLL_EXPORT void dsp_file_write_png_composite(const char* filename, int components, int compression, dsp_stream_p* stream);
 
DLL_EXPORT dsp_t* dsp_file_bayer_2_gray(dsp_t* src, int width, int height);
 
DLL_EXPORT dsp_t* dsp_file_bayer_2_rgb(dsp_t *src, int red, int width, int height);
 
DLL_EXPORT dsp_stream_p *dsp_stream_from_components(dsp_t* buf, int dims, int *sizes, int components);
 
DLL_EXPORT dsp_stream_p *dsp_buffer_rgb_to_components(void* buf, int dims, int *sizes, int components, int bpp, int stretch);
 
DLL_EXPORT void dsp_buffer_components_to_rgb(dsp_stream_p *stream, void* rgb, int components, int bpp);
 
DLL_EXPORT dsp_t* dsp_file_composite_2_bayer(dsp_stream_p *src, int red, int width, int height);
 
DLL_EXPORT void dsp_file_write_fits_bayer(const char* filename, int components, int bpp, dsp_stream_p* stream);
 
DLL_EXPORT dsp_t* dsp_file_bayer_2_composite(dsp_t *src, int red, int width, int height);
 
DLL_EXPORT dsp_align_info *dsp_align_fill_info(dsp_triangle t1, dsp_triangle t2);
 
DLL_EXPORT dsp_triangle *dsp_align_calc_triangle(dsp_star* stars, int num_stars);
 
DLL_EXPORT void dsp_align_free_triangle(dsp_triangle *triangle);
 
DLL_EXPORT int dsp_align_get_offset(dsp_stream_p ref, dsp_stream_p to_align, double tolerance, double target_score, int num_stars);
 
#include <fits_extensions.h>
 
#ifdef  __cplusplus
}
#endif
 
#endif //_DSP_H