Вопрос по полиномиальной интерполяции

[Бахурин Сергей]

В принципе да. Но проблемы стыковки будут те же при блочной обработке. При этом поднять частоту в 160 раз видится довольно сложной. Лучше один раз отладить функцию блочного ресамплинга

[Бахурин Сергей]

Вот программа, которая делает пересчет и умеет правильно стыковать блоки

Код: Выделить всё

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"

// длина исходного сигнала
#define N   256

// размер блока
#define L   16

// P/Q
#define P   160
#define Q   147


int main(int argc, char* argv[])
{
    void* hdspl;  /* DSPL handle        */
    void* hplot;  /* GNUPLOT handle     */

    double s[N], ts[N];   // исходный сигнал на входе ресамплера и время
    double dt;            // интервал дискретизации после ресамплера Q/P/Fs
    double *ty = NULL;    // Время после ресамплинга
    double *y = NULL;     // Ресамплинг в блочном режиме длины L = 16 
    double *z = NULL;     // Ресамплинг для всего сигнала

    int k, ny;
    int err;

    /* Load DSPL function  */
    hdspl = dspl_load();

    // Исходный сигнал
    for(k = 0; k < N; k++)
    {
        ts[k] = (double)k;
        s[k] = sin(M_2PI*0.2176870748*ts[k]);
    }
    // интервал дискретизации после ресамплера
    dt = (double)Q/(double)P;
    
    // Количество отсчетов после ресамплера
    ny = (int)((double)(N-1)/dt)+1;
    
    y  = (double*) malloc(ny * sizeof(double));
    ty = (double*) malloc(ny * sizeof(double));
    
    // Время после ресамплера (для вывода на график)
    for(k = 0; k < ny; k++)
        ty[k] = (double)k*dt;
    
    // преременные пересчета времени в массивах s и y
    // frd - временное смещение отсчета на выходе ресамплера отностиельно входа
    double ts0,ts1, ty0, ty1, frd;
    int pos, ypos; // позиции в массивах
    ts0 = 0.0;
    ty0 = 0.0;
    ts1 = ty1 = 0.0;
    
    frd = 0.0;
    pos = 0;
    ypos = 0;
    
    // выход блока
    double* tmp = NULL;
    int ntmp;
    
    while(pos+L < N)
    {
        // обрабатываю текущий блок начиная с позиции pos (L - Длина блока)
        farrow_spline(s+pos, L, P, Q, frd, &tmp, &ntmp);
        
        // копирую в выходной массив со сдигом 2 отсчета, 
        // чтобы убрать эффекты стыковки блоков
        if(ypos)
            memcpy(y+ypos+2, tmp+2, (ntmp-2)*sizeof(double));
        else
            memcpy(y, tmp, (ntmp)*sizeof(double));

        // Пересчитываю времена конца блоков для входа и выхода
        ty1 = ty0+(double)ntmp*dt;
        ts1 = ts0+(double)L;
        printf("pos = %d, ts1 = %.3f, ty1 = %.3f\n", pos,ts1,ty1);
        
        // позиции на выходе ресамплера на след блок
        ypos += ntmp;
        // делаю  смещение на 4 отсчета раньше, чтобы было перекрытие между блоками
        ts0  = ts1 - 4.0;
        
        // теперь надо найди начало массива после ресамплинга, 
        // чтобы посчитать сдвиг во времени
        ty0 = ty1;
        while(ty0>ts0)
        {
            ty0-=dt;
            ypos--;
        }
        
        // считаю сдвиг во времни от 0 до 1
        frd = ts0 - ty0; 
        printf("ts0 = %.3f, ty0 = %.3f frd = %.3f\n", ts0,ty0, frd);
        
        // Смещение позиции на 4 отсчета раньше для перекрытия
        pos += L-4;
    }
    farrow_spline(s+pos, N-pos, P, Q, frd, &tmp, &ntmp);
    memcpy(y+ypos+2, tmp+2, (ntmp-2)*sizeof(double));

    // Ресамплинг для всего куска сразу
    farrow_spline(s, N, P, Q, 0, &z, &k);

    /* save input signal and filter output to the txt-files */
    // если открыть файлы, будет видно, что массивы y и z полностью совпадают
    writetxt(ts, s, N,  "dat/s.txt");
    writetxt(ty, y, ny, "dat/y.txt");
    writetxt(ty, z, ny, "dat/z.txt");

    /* plotting by GNUPLOT*/
    gnuplot_create(argc, argv, 820, 540, "img/resampling_block.png", &hplot);
    gnuplot_cmd(hplot, "unset key");
    gnuplot_cmd(hplot, "set grid");
    gnuplot_cmd(hplot, "set xlabel 'n'");
    gnuplot_cmd(hplot, "set ylabel 's(n)'");
    gnuplot_cmd(hplot, "set yrange [-1.5:1.5]");
    gnuplot_cmd(hplot, "set multiplot layout 2,1 rowsfirst");
    gnuplot_cmd(hplot, "plot 'dat/s.txt'  with lines");

    gnuplot_cmd(hplot, "plot 'dat/y.txt' with lines,\\");
    gnuplot_cmd(hplot, "     'dat/z.txt' with lines");
    gnuplot_cmd(hplot, "unset multiplot");
    gnuplot_close(hplot);
    
    /* free DSPL handle */
    dspl_free(hdspl);

    if(ty)
        free(ty);
    if(y)
        free(y);
    if(z)
        free(z);
    if(tmp)
        free(tmp);
    return err;
}

Если выполнить на верхнем графике исходный сигнал, на нижнем ресамплинг блочный и для всего исходного массива совпадают.

[alexx_78]

Понял, спасибо! Буду изучать