передискретизация, как? Опции

[TigerSHARC]

Здравствуйте.
Есть сигнал:

Fs = 6400;
N = 128;
Ts = 1/Fs;
T = Ts*N;
f = 45;
t = 0:Ts:T-Ts;
A = 10;
Y = A*sin(2*pi*f*t)+ A*sin(2*pi*2*f*t)+ A*sin(2*pi*3*f*t);

как и какой функцией в MATLAB пересчитать частоту дискретизации с 6400 на 5760.
Причём с тем же числом отсчётов (128).

Это нужно для корректировки спектра сигнала.
Дело в том, что когда меняю частоту "руками":

Fs = 5760;
N = 128;
Ts = 1/Fs;
T = Ts*N;
f = 45;
t = 0:Ts:T-Ts;
A = 10;
df = Fs/N;
f1 = 0:df:Fs-(Fs/N);
Y = A*sin(2*pi*f*t)+ A*sin(2*pi*2*f*t)+ A*sin(2*pi*3*f*t);
X = fft(Y)/length(t);
>> stem(f1, abs(X)), grid

получается красивый спектр, отражающий верную информацию об амплитудах сигнала.
Пробовал для передискреизации функцию resample, но она выдаёт меньшее число отсчётов и спектр искажён...

да и вообще, не толкьо в MATLAB, как из 128 отсчётов получается 128 передискретизированных отсчётов... полюбому время выборки должно как-то варьироваться...

[TigerSHARC]

тогда так.
1) при интерполяции всегда используется фиксированное число отсчётов с АЦП?
2) количесвто отсчётов(идущих на блок БПФ) после ресемплинга остаётся прежним? (в моём случае 128)

Просто использую я и фильтр фарроу на С и матлаб... ну не получается и всё... когда руками меняю частоту

Fs = 5760;
N = 128;
Ts = 1/Fs;
T = Ts*N;
f = 45;
t = 0:Ts:T-Ts;
A = 10;
df = Fs/N;
f1 = 0:df:Fs-(Fs/N);
Y = A*sin(2*pi*f*t)+ A*sin(2*pi*2*f*t)+ A*sin(2*pi*3*f*t);
X = fft(Y)/length(t);
>> stem(f1, abs(X)), grid

тогда всё окей.
но тогда и меняется и длинна окна (стала 0,0222).
как это при ресемплинге происходит?

[bahurin]

как это при ресемплинге происходит?

Если частота понижается, то для сохранения количества отсчетов ПОСЛЕ ресамплинга надо уведить ИСХОДНОЕ количество отсчетов, потом применить рэсамплер.

Код

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define M_PI       3.14159265358979323846

// Функция сохранения результатов анализа в текстовый файл для
// построения графиков.
void saveToTXT(double* t, double* s, int n, char* fileName){
    FILE *file = fopen(fileName, "w");
    if(!file) return;
    for(int i = 0; i < n; i++)
        fprintf(file,"%e\t%e\n",t[i],s[i]);
    fclose(file);
}


// Функция расчета кубического полинома
// на основе модифицированного фильтра Фарроу
double farrow3(double *y, double x){
    double a3 = (y[3]-y[0])/6.0+(y[1]-y[2])/2.0;
    double a1 = (y[3]-y[1])/2.0-a3;
    double a2 = y[3]-a3-a1-y[2];
    return x*(x*(x*a3+a2)+a1)+y[2];
}

int main(){
            
    double Fs0 = 16E3;            // исходная частота дискретизации
    double Fs1 = 12E3;            // частота дискретизации после ресэмплинга
    
    int K = 128;// количество отсчетов после ресэмплинга
    int N = (int)K*Fs0/Fs1;// количество исходных отсчетов
    // выделяю память под исходный сигнал
    double* y0 = (double*) malloc(N*sizeof(double));
    double* t0 = (double*) malloc(N*sizeof(double));
    // выделяю память под результат ресэмплинга
    double* y1 = (double*) malloc(K*sizeof(double));
    double* t1 = (double*) malloc(K*sizeof(double));
    // формирую исходный сигнал
    for(int i = 0; i<N; i++){
        t0[i] = double(i)/Fs0;
        y0[i] = cos(2*M_PI*3E3*t0[i]);        // гармоническое колебание на частоте 3 кГц.
    }
    // произвожу ресэмплинг
    for(int i = 0; i<K; i++){
        t1[i] = (double)i/Fs1;                // текущий момент дискретизации
        int j = int((double)i*Fs0/Fs1)-2;    // пересчитываю индекс исходного сигнала
        if(j<0) j =0;                        // если индекс отрицателььный то приравниваю его к 0
        double x = (t1[i]-t0[j])*Fs0-2;        // пересчитываю значение x
        y1[i] = farrow3(y0+j, x);            // фильтр Фарроу
        t1[i] *= 1E3;                        // перевожу время результата ресэмплинга в мс
    }
    // перевожу исходное время в мс
    for(int i = 0; i<N; i++){
        t0[i] *= 1E3;                        
    }
    // сохраняю результат в файлы
    saveToTXT(t0, y0, N, "исходный сигнал.txt");
    saveToTXT(t1, y1, K, "resampling.txt");
    system("Pause");
    return 0;
}

все работает.

Сообщение отредактировал bahurin - Feb 22 2010, 17:40

[TigerSHARC]

Если частота понижается, то для сохранения количества отсчетов ПОСЛЕ ресамплинга надо уведить ИСХОДНОЕ количество отсчетов, потом применить рэсамплер.

Код

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define M_PI       3.14159265358979323846

// Функция сохранения результатов анализа в текстовый файл для
// построения графиков.
void saveToTXT(double* t, double* s, int n, char* fileName){
    FILE *file = fopen(fileName, "w");
    if(!file) return;
    for(int i = 0; i < n; i++)
        fprintf(file,"%e\t%e\n",t[i],s[i]);
    fclose(file);
}


// Функция расчета кубического полинома
// на основе модифицированного фильтра Фарроу
double farrow3(double *y, double x){
    double a3 = (y[3]-y[0])/6.0+(y[1]-y[2])/2.0;
    double a1 = (y[3]-y[1])/2.0-a3;
    double a2 = y[3]-a3-a1-y[2];
    return x*(x*(x*a3+a2)+a1)+y[2];
}

int main(){
            
    double Fs0 = 16E3;            // исходная частота дискретизации
    double Fs1 = 12E3;            // частота дискретизации после ресэмплинга
    
    int K = 128;// количество отсчетов после ресэмплинга
    int N = (int)K*Fs0/Fs1;// количество исходных отсчетов
    // выделяю память под исходный сигнал
    double* y0 = (double*) malloc(N*sizeof(double));
    double* t0 = (double*) malloc(N*sizeof(double));
    // выделяю память под результат ресэмплинга
    double* y1 = (double*) malloc(K*sizeof(double));
    double* t1 = (double*) malloc(K*sizeof(double));
    // формирую исходный сигнал
    for(int i = 0; i<N; i++){
        t0[i] = double(i)/Fs0;
        y0[i] = cos(2*M_PI*3E3*t0[i]);        // гармоническое колебание на частоте 3 кГц.
    }
    // произвожу ресэмплинг
    for(int i = 0; i<K; i++){
        t1[i] = (double)i/Fs1;                // текущий момент дискретизации
        int j = int((double)i*Fs0/Fs1)-2;    // пересчитываю индекс исходного сигнала
        if(j<0) j =0;                        // если индекс отрицателььный то приравниваю его к 0
        double x = (t1[i]-t0[j])*Fs0-2;        // пересчитываю значение x
        y1[i] = farrow3(y0+j, x);            // фильтр Фарроу
        t1[i] *= 1E3;                        // перевожу время результата ресэмплинга в мс
    }
    // перевожу исходное время в мс
    for(int i = 0; i<N; i++){
        t0[i] *= 1E3;                        
    }
    // сохраняю результат в файлы
    saveToTXT(t0, y0, N, "исходный сигнал.txt");
    saveToTXT(t1, y1, K, "resampling.txt");
    system("Pause");
    return 0;
}

все работает.

у них там да. А вы задайте частоту дискретизации 6400 для синуса c частотой f0 = 45,35Гц. при выборке N = 1024 (такое соотношение частоты дискретизации и количества выборок дают верный спектр лишь при f0 = 50Гц, так как тогда в окне ровно 8 периодов), а мы же для f0 = 45,35Гц здесь видим нецелое число периодов в окне.
затем попробуйте передискретизировать этот синус на на частоту 5804.8 (заметьте дробная частота получается как N/8*f0) и должны увидеть целое количесвто периодов в окне.

Но в моих экспериментах с этим кодом получается ерунда при передискретизации.

Ну бог с ним, с эти кодом С-шным. Хотябы в MATLAB бы увидеть положительный результат.

[bahurin]

Но в моих экспериментах с этим кодом получается ерунда при передискретизации.

ставьте

Код

    

double Fs0 = 6400;            // исходная частота дискретизации
double Fs1 = 5804.8;        // частота дискретизации после ресэмплинга
    
int K = 1024;// количество отсчетов после ресэмплинга
int N = (int)K*Fs0/Fs1;// количество исходных отсчетов

а также сигнал:

Код

    

// формирую исходный сигнал
    for(int i = 0; i<N; i++){
        t0[i] = double(i)/Fs0;
        y0[i] = cos(2*M_PI*45.35*t0[i]);        // гармоническое колебание на частоте 3 кГц.
    }

все замечательно работает. Количество входных отсчетов 1128 выходный 1024. На выходе ровно 8 периодов колебания