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

[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 передискретизированных отсчётов... полюбому время выборки должно как-то варьироваться...

[Andron_]

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

а как это с точки зрения такой науки как арифметика возможно?
длительность сигнала одна и та же, число отсчетов в единицу времени меньше. откуда дополнительным отсчетам взяться?

что вы хотите добиться - я не понял.

[bahurin]

Передискретизацию можно выполнить на основе интерполятора. Часто прибегают к лагранжевой интерполяции и фильтрам фарроу. Там же есть программная реализация на с++. При ресамплигнге на более низкую частоту дискретизации не забудьте поставить фильтр антиалисинга. Однако если вы хотите получить спектр сигнала, то вам стоит рассмотреть окна сглаживания

[TigerSHARC]

с точки зрения такой науки как численные методы возможно провести интерполяцию сигнала.

P.S. не использую оконное взвешиване здесь так как ГОСТ это запрещает (для моего случая) вот и приходится частоту менять программно.

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

[bahurin]

P.S. не использую оконное взвешиване здесь так как ГОСТ это запрещает (для моего случая)

В штате Индиана (США) в 1897 году был выпущен билль (см.: en:Indiana Pi Bill), законодательно устанавливающий значение числа Пи равным 3,2.

Вобще странная ситуация, когда гост запрещает чтото использовать. Почему бы не использовать окна, если они позволяют получить установленные гостом харакетристики? Чудно однако.

[Andron_]

чтобы число отсчетов осталось тем же, нужна экстраполяция.

[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]

В штате Индиана (США) в 1897 году был выпущен билль (см.: en:Indiana Pi Bill), законодательно устанавливающий значение числа Пи равным 3,2.

Вобще странная ситуация, когда гост запрещает чтото использовать. Почему бы не использовать окна, если они позволяют получить установленные гостом харакетристики? Чудно однако.

не хотелось бы чтобы дисскусия ушла в рксло обсуждения гостов, НО...!

Там так сказано видимо потому, что частотное разрешение должно быть 5 Гц ( а не 50 как раньше)*, чтобы интергармоники (кратные 5 Гц) содержали полезную информацию, а не информацию о расползаном спектре.
А чтобы обеспечить такое разрешние по частоте придётся брать выборку с огромным запасом, что скажется на быстродействии в реальном времени (к тому же предполагается работа сразу 8-ми каналов одновременно по 3 фазы тока и напряжения и нейтраль)
Алгоритм с использованием окна, по-моему, на 8 каналов DSP проц не потянет (по моему)...
Подход с передискретизацией и прямоуг. окне хорош тем, что позволит при заданной точности обеспечить минимальную выборку сигнала.
Я думаю так.

* имеется ввиду гост на качество электроэнергии и нововведённые госты на измерение гармоник и интергармоник в электрических сетях.



но я не о ГОСТах я о самом алгоритме. Возможно же такое преобразование сигнала, чтобы в окно уместилось целое (с некоторой малой погрешностью) количесвто периодов. НО как!?

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

... иными словами я хочу "поджать сигнал" во временной области под окно 0,02с.(даже если его период больше, или меньше этого значения), сохраненив его параметы. При передискретизации происходит этот эффект?



но я не о ГОСТах я о самом алгоритме. Возможно же такое преобразование сигнала, чтобы в окно уместилось целое (с некоторой малой погрешностью) количесвто периодов. НО как!?

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

... иными словами я хочу "поджать сигнал" во временной области под окно 0,02с.(даже если его период больше, или меньше этого значения), сохраненив его параметы. При передискретизации происходит этот эффект?



но я не о ГОСТах я о самом алгоритме. Возможно же такое преобразование сигнала, чтобы в окно уместилось целое (с некоторой малой погрешностью) количесвто периодов. НО как!?

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

... иными словами я хочу "поджать сигнал" во временной области под окно 0,02с.(даже если его период больше, или меньше этого значения), сохраненив его параметы. При передискретизации происходит этот эффект?

[bahurin]

не хотелось бы чтобы дисскусия ушла в рксло обсуждения гостов, НО...!

Согласен воздух сотрясать глупо. Выше программа на си которая делает ресамплинг с фиксированным количеством отсчетов на выходе

[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 периодов колебания

[TigerSHARC]

ставьте

Код

    

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 периодов колебания

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

Да и почему то исходное количесвто отсчётов разное всегда... наверно исходное число должно быть постоянным а выходное менятся в зависимости от соотношения F0/F1. Так?

Да и скажите где вы этот код взяли.
тот что у меня с того же сайта:

int main()
{
int N = 64; // количество исходных отсчетов
double Fs0 = 20E3; // исходная частота дискретизации
double Fs1 = 15E3; // частота дискретизации после ресэмплинга
int K = (int)(N*Fs1/Fs0); // количество отсчетов после ресэмплинга
// выделяю память под исходный сигнал
.......


вот тут как раз исходное количесвто отсчётов постоянно(как и должно быть), а на выходе меняется...

странно

Сообщение отредактировал TigerSHARC - Feb 23 2010, 13:32

[Andron_]

это что-то нереальное...

сдайте дилера, я тоже хочу...

[TigerSHARC]

Здаётся мне у меня какие-то принципиальные затруднения....
Хочу просто "поджать" сигнал во временной области под жестко заданное временное окно с сохранением информации о сигнале.
Гуру, подскажите!!!!

[Andron_]

вы не можете его "поджать" под временное окно. Не, при желании, конечно, можно, но при этом спектр изменится. Сигнал можно только _обрезать_ под временное окно. О целом числе периодов сигнала на данном временном окне речи не идет.
Можно подогнать сигнал определенной длительности под нужное число отсчетов изменением частоты дискретизации.

[TigerSHARC]

вы не можете его "поджать" под временное окно. Не, при желании, конечно, можно, но при этом спектр изменится. Сигнал можно только _обрезать_ под временное окно. О целом числе периодов сигнала на данном временном окне речи не идет.
Можно подогнать сигнал определенной длительности под нужное число отсчетов изменением частоты дискретизации.

Ну так я и хочу всего лишь навсего чтобы частота чигнала точно легла на бины Фурье. хочу добиться того, чтобы соотношение частоты дискретизации и длинны выборки - (т.е. частотное разрешение фурье), позволило увидеть спектр сигнала в виде одной палки, а не гребёнки от растекания.

P.S. Окна не предлагать. Рочему? См. выше.

[Andron_]

тью... "вон оно че, Михалыч"...

тогда, как я понимаю, 2 этапа.
1. выбрать длительность сигнала, на которую уложится целое число периодов сигнала.
2. выбрать новую частоту дискретизации таковой, чтобы на выбранное время приходилось 2^n отсчетов.


а, ну и третий - в матлабе применить resample на подсчитанную частоту дискретизации.

[TigerSHARC]

тью... "вон оно че, Михалыч"...

тогда, как я понимаю, 2 этапа.
1. выбрать длительность сигнала, на которую уложится целое число периодов сигнала.
2. выбрать новую частоту дискретизации таковой, чтобы на выбранное время приходилось 2^n отсчетов.


а, ну и третий - в матлабе применить resample на подсчитанную частоту дискретизации.

resample не канает(((

пример

Fs = 6400;
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);
Yr = resample(Y, 9,10);
subplot(2,1,1), stem (t,Y), grid
subplot(2,1,2), stem (t, Yr), grid

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

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, но она выдаёт меньшее число отсчётов и спектр искажён...


функцию resample, но она выдаёт меньшее число отсчётов и спектр искажён...

Сообщение отредактировал TigerSHARC - Feb 23 2010, 19:40

[leksa]

Здравствуйте!
Может быть вот такой вариант Вас устроит:
% variant
% накапливаем сигнал на исходной частоте
Fs = 6400;
N = 143;
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);
% интерполяция и БПФ
N2 = 128
Fs2 = 5760
Ts2 = 1/Fs2
ti = [0:Ts2:(N2-1)*Ts2]
Yi = interp1(t,Y,ti,'spline');% интерполяция методом сплайнов, можно и другие методы, например, полиномиальную интерполяцию 2 и 3 порядка

figure
plot(t,Y,'b- x'),grid on,hold on
plot(ti,Yi,'r- o')
Xi = fft(Yi)/length(ti);
df2 = Fs2/N2;
f2 = 0:df2:Fs2-(Fs2/N2);
figure
stem(f1, abs(X)), grid on,hold on
figure
stem(f2, abs(Xi),'g o'), grid

%>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

Если вы хотите изменить частоту дискретизации с 6400 на 5760 так, чтобы на частоте 5760 у Вас оказалось 128 отсчетов, то для интерполяции с частоты 6400 вам понадобится кусок сигнала длительностью 128*6400/5760 ~= 143 отсчета.
Взять 128 отсчетов на частоте 6400 и переделать их в 128 отсчетов на частоте 5760 невозможно, так как длительность вашего сигнала на частоте 6400 равна 128/6400 = 0.0200 сек, а длительность сигнала в 128 отсчетов на частоте 5760 равна 0.0222 секунды.
Где взять недостающих отсчетов сигнала на 0.0022 секунды?
2 варианта:
1 - накопить на частоте 6400 (скрипт выше в этом посте).
2 - заниматься экстраполяцией, то есть пытаться продолжить сигнал в будущее на основе имеющихся отсчетов. Успех такой операции зависит от метода экстраполяции и самого экстраполируемого сигнала. Первый способ я думаю проще и надежнее.

Сообщение отредактировал leksa - Feb 23 2010, 21:09

[TigerSHARC]

на самом деле нужно просто передискретизировать сигнал так, чтобы частота дискретизации была кратна частоте синуса. И при этом частота точно ляжет на бины фурье. вот и всё.
Я знаю, что для этих целей применяется фильтр Фарроу или полифазные КИХ фильтры (или просто фильтры дробной задержки).
На С пробовал фильтр фарроу - результата нет. Спектр тот же примерно получается после передискретизации.
А нужно так - получили синус число периодов которого не уместилось во временном окне, а спектр вместо одной палки на частоте синуса содержит гребёнку на всех частотах. Затем произвели передискретизацию. И вот число периодов целое во временном окне, а спектр уже одна палка на частоте синуса, а его амплитуда известна с какой-то малой погрешностью...

leksa, в любом случае огромное спасибо за код!

Сообщение отредактировал TigerSHARC - Feb 23 2010, 21:30

[leksa]

на самом деле нужно просто передискретизировать сигнал так, чтобы частота дискретизации была кратна частоте синуса. И при этом частота точно ляжет на бины фурье. вот и всё.
Я знаю, что для этих целей применяется фильтр Фарроу или полифазные КИХ фильтры (или просто фильтры дробной задержки).
На С пробовал фильтр фарроу - результата нет. Спектр тот же примерно получается после передискретизации.
А нужно так - получили синус число периодов которого не уместилось во временном окне, а спектр вместо одной палки на частоте синуса содержит гребёнку на всех частотах. Затем произвели передискретизацию. И вот число периодов целое во временном окне, а спектр уже одна палка на частоте синуса, а его амплитуда известна с какой-то малой погрешностью...

Я понял что вы это хотите сделать.
Приведенный мной выше код позволяет это сделать, с одним условием накапливать на более высокой частоте отсчетов нужно больше, чем их должно быть на низкой. Какой конкретно интерполятор использовать - не так важно, имхо. Но полифазные КИХ работают только если изменение частоты выражается соотношением целых чисел, 5760/6400=47/80 - несократимая дробь, имхо, при таком соотношении полифазная структура не очень эффективна, хотя точно сейчас не скажу.

[leksa]

За код - пожалуйста, тем более это ваш код в общем то с небольшими изменениями.
По полифазным КИХ я выше написал неверно - вполне их тоже в этой ситуации можно использовать.
Извините, 5760/6400 = 9/10 конечно. В матлабе факторизовал, забил 3760 вместо 5760 - по ночам наверно лучше спать! )

вот вариант с полифазными КИХ
% variant
Fs = 6400;
N = 142;
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);

N2 = 128
Fs2 = 5760
Ts2 = 1/Fs2
ti = [0:Ts2:(N2-1)*Ts2]
Yi = resample(Y,9,10);%interp1(t,Y,ti,'spline');

figure
plot(t,Y,'b- x'),grid on,hold on
plot(ti,Yi,'r- o')
Xi = fft(Yi)/length(ti);
df2 = Fs2/N2;
f2 = 0:df2:Fs2-(Fs2/N2);
figure
stem(f1, abs(X)), grid on,hold on
figure
stem(f2, abs(Xi),'g o'), grid

Сообщение отредактировал leksa - Feb 23 2010, 22:53

[Andron_]

resample не канает(((

пример

Fs = 6400;
N = 128;
...

вы иногда думайте над тем, что вам советуют.

КАК вы собрались из 128-ми отсчетов исходного сигнала на более высокой частоте дискретизации получить столько же, но на более низкой?
Законы мироздания таковы, что это слабо выполнимая задача.

Считайте.
128 отсчетов. частота дискретизации 5760. Длительность сигнала 128/5760 = 0,022(2)
Число отсчетов, необходимое для такой длительности сигнала при частоте дискретизации 6400: 6400 * 0,022 = 142,22(2) ~ 143 отсчета.

Соответственно и в коде своем пишите:

Fs = 6400;
N = 143; %!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
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);
Yr = resample(Y, 9,10);
subplot(2,1,1), stem (t,Y), grid
subplot(2,1,2), stem (t, Yr), grid

[Andron_]

только в коде нужно сначала передискретизацию сделать, а затем фурье.

[TigerSHARC]

дело в том что нужно обеспечить отсутсвтие растекания спектра, если основная гармоника колеблется в пределах 42.5...57.5Гц
теперь ясно делаю так:
послольку самая низкая частота 42.5, то нужно всегда брать временное окно равное одному периоду этой частоты, так как все остальные периоды в диапазоне частот до 57.2 укладываются целиком на этом интервале.
Затем беру 128 точек на такой частоте дискретизации 128*f0. Получается как бы вырезаем ровно период из сигнала.
Только придёться алгоритм определения частоты искать (что выходит за рамки данной темы).
Априори частота полученя с абсолютной точностью.
Вот код:

Fs = 6400;
Ts = 1/Fs;
N = (1/42.5)/Ts;
N = round((1/42.5)/Ts);
T = Ts*N;
f = 55;
t = 0:Ts:1/42.5-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 = abs(fft(Y)/length(t));
% интерполяция и БПФ
N2 = 128;
Fs2 = 128*f; % этот параметр выбирается исходя из найденной по частотному алгоритму частоты основной гармоники
Ts2 = 1/Fs2;
ti = [0:Ts2:(N2-1)*Ts2];
Yi = interp1(t,Y,ti,'spline');
Xi = abs(fft(Yi)/length(Yi));

сравниваем спектры - всё классно!

Жду коментарии.

Сообщение отредактировал TigerSHARC - Feb 24 2010, 10:36

[TigerSHARC]

Или получилось из пустого в порожнее? по сути просто вырезал кусок сигнала.
что не желательно система должна оценивать сигнал непрерывно.
так вот, к сожалению буду голословен, но должно быть примерно так:
...интерполятор работает так что берет такое количество отсчетов что в них содержится ровно один период, независимо от частоты.

И это действительно так. Реализуется с помощью фарроу-фильтра.

незнаю как пока что...

Сообщение отредактировал TigerSHARC - Feb 24 2010, 12:33

[SPACUM]

Реализуется с помощью фарроу-фильтра.

[/quote]
А сколько гармоник надо? И какая точность? Если всего 3, тогда DFT любой FFT обгонит в несколько раз даже без Фарроу. И в несколько раз по точности. И никакого ресемплинга не требуется. А если много и уровни маленькие, тогда Фарроу там такое наделает... Попробуйте подставить реальные уровни гармоник, погрешность у Фарроу может ползти от уровня основной гармоники.
Я делал так:
1. Начало измерений от прохождения сигнала через 0.
2. Конец измерений до прохождения сигнала через 0 после требуемого числа периодов.
3. Произвести требуемые DFT на столько точек, сколько попадет в период измерения.
Какое-нибудь окно, Ханна например, повысит точность. Измерять частоту не надо. Можно в реале через период, если проц быстрый.
Попробуйте, гармоники должны быть реальные не более 5-10%.
Интересно какой ГОСТ допускает ресэмплинг?

[TigerSHARC]

Да, то что у меня получилось выше просто вырезает период сигнала. Интрвал неопределённости в моём случае (если основная гармоник пляшет 42,5....57,5Гц) составит порядка 0,006 с. (разность периодов)
Если применить окно ханна получается супер точно... но его прменять нельзя!!!!!
Давайте расскажу всё как есть, дабы прояснить ситуацию.
Создаётся алгоритм по ГОСТ МЭК 61000-4-7.
Аппаратно частоту дискретизации менять нельзя.
Ресемплинг применяется уже в разработанном устройстве вот ссылка: http://powerdsp.narod.ru/analizatorkachestva/measpar.pdf

если не лень прочитайте.
Там применяется большая частота дискретизации и выборка длиннее.
Цитирую разработчика:
"... количество отсчетов фурье фиксировано но на него приходится строго один период, интерполятор работает так что берет такое количество отсчетов что в них содержится ровно один период, независимо от частоты."

Общался с разработчиком - для 6400Гц и 128 выборок ресемплинг тоже вполне уместен.

Я пишу алгоритм на С. Получается тольуо "вырезать" сигнал как в коде выше...(((

Сообщение отредактировал TigerSHARC - Feb 24 2010, 15:13

[leksa]

Я прочитал документ по ссылке.
Из того что вы писали все равно не ясно, чего вы хотите добиться.
В приведенном вами документе сначала производится оценка частоты первой гармоники сигнала путем подсчета переходов через ноль фильтрованного режекторным фильтром сигнала. Затем используя эту оценку, интерполятор выдает отсчеты сигнала так, чтобы на размер ДПФ пришлось целое число периодов, тогда все гармоники тоже "лягут" целое число раз.
Делайте также. Или что вас не устраивает?

[TigerSHARC]

Уважаемый, leksa, если вы поняли суть алгоритма - поделитесь пожалуйста соображениями.

Я именно это и пытаюсь сделать. Только у меня, в приведённом выше коде, происходит как бы вырезание одного периода из временного окна. В результате во временном окне остаётся "не тронутой" небольшая область. Эта область тем больше, чем больше частота сигнала (максимальная когда f = 57.5).
А это недопустимо, так как по госту паузы между временными окнами недопускаются....

Вобщем тробный ресемплер (в частности фарроу) должен привести к тому, что окно наблюдения будет менять свою длительность в зависимости от частоты основной гармоники.

[TigerSHARC]

Уважаемый, leksa, если вы поняли суть алгоритма - поделитесь пожалуйста соображениями.

Я именно это и пытаюсь сделать. Только у меня, в приведённом выше коде, происходит как бы вырезание одного периода из временного окна. В результате во временном окне остаётся "не тронутой" небольшая область. Эта область тем больше, чем больше частота сигнала (максимальная когда f = 57.5).
А это недопустимо, так как по госту паузы между временными окнами недопускаются....

Вобщем тробный ресемплер (в частности фарроу) должен привести к тому, что окно наблюдения будет менять свою длительность в зависимости от частоты основной гармоники.

Сообщение отредактировал TigerSHARC - Feb 24 2010, 16:58

[thermit]

Дык, пауз во времени и не будет. Меняться будет число отсчетов в окне.

[TigerSHARC]

Дык, пауз во времени и не будет. Меняться будет число отсчетов в окне.

уменя выходит, что например в исходном сигнале за периодом основной частоты наблюдается сильный всплеск. если в овно умещается ровно период и ещё кусок этого всплеска, то после интерполяции получаю просто период, а о всплеске за периодом информация теряется...
И число отсчётов в окне изначально у меня фиксированно а после ресемплинга всегда равно 128 (но соответсвенно для разной частоты дискретизации).... Вобщем получается что после ресемплинга вырезаю просто период.
Довольно трудно объяснить такие вещи, но думаю у меня получилось)

расскажите в чём неправ...

[TigerSHARC]

Вобщем, прошу подсказать как сделать адаптивное временное окно на основе передискретизации сигнала.
(формулировка по-моему правильна)
Причём берём каждый раз постоянное число отсчётов с АЦП.

Сообщение отредактировал TigerSHARC - Feb 24 2010, 18:24

[leksa]

TigerSHARC
Извините, но мне все-таки не до конца ясны все условия вашей задачи, да и потом предлагать какие-то конкретные решения - их нужно проверять на практике применительно к вашей задаче, то есть сделать, а это время.
Но все же, вот какие моменты мне до сих пор не ясны -
1) Какой интервал наблюдения, за который вы должны получить оценку уровня сигнала?
2) Требуется оценить уровень первой гармоники, или еще какие-то параметры?
3) Какой период обновления полученных оценок? То есть через какой интервал времени нужно обновлять полученные оценки?
Если интервал наблюдения - это один период первой гармоники, то может быть делать так:
Предположим на вход алгоритма поступают отсчеты с частотой дискретизации 6400 Гц.
Алгоритм отслеживает переходы во входном сигнале через ноль, причем при переходе от отрицательного уровня к положительному, запускается накопление сигнала в некий буфер.
При следующем переходе через ноль от отрицательного уровня к положительному, накопление в буфер заканчивается.
На этот момент у вас в буфере 1 период синуса.
По числу накопленных отсчетов можно оценить частоту первой гармоники, правда довольно грубо так как погрешность оценки периода у вас +/- 1 отсчет. Может быть, оценку можно улучшить, если момент перехода через ноль уточнять линейной интерполяцией.
То есть у вас один отсчет был ниже нуля, следующий выше. Так как их уровни вам известны, то можно построить уравнение прямой линии через эти два отсчета и определить в какой момент линия пересекает нулевой уровень. Так нужно делать два раза в начале и конце периода. Потом период вычислите как разницу моментов перехода через ноль.
Теперь период сигнала известен. Вычисляем массив времени t=0:T1/N:(T1/N)*(N-1), где T1 - оценка периода первой гармоники, а N - размерность ДПФ, 128 например.
Используя этот массив времени и накопленный в буфере сигнал интерполятором получаем те N отсчетов одного периода, которые вам нужны для вычисления ДПФ, чтобы на интервал ДПФ попал точно 1 период.
По уровню первой гармоники ДПФ получаем искомую оценку уровня сигнала с частотой первой гармоники.
Сбрасываем содержание буфера и начинаем накопление заново, до нового перехода через ноль сигнала по направлению от отрицательных значений к положительным.
В таком варианте, оценка будет обновляться периодически, но не через 128 отсчетов входного сигнала, а через 1 период первой гармоники. Будет гармоника 45 Гц - будет реже обновляться, 52 Гц -чаще.
Данная схема чисто из головы, так что я не претендую на то что это наилучший способ.
Чтобы такое в матлабе набросать требуется время, а его мало.

Да, и при высоком уровне помех/искажений, моменты перехода через ноль будут "дрожать" слишком сильно, то есть оценка периода будет не точной, тогда и оценка уровня первой гармоники будет искажена. Как быстро может меняться частота первой гармоники? Если скажем на доли Гц в минуту, то тогда может быть все-таки накапливать побольше периодов, по переходам через ноль определять частоту первой гармоники. Такая оценка будет менее подвержена влиянию дрожания переходов через ноль.

Сообщение отредактировал leksa - Feb 24 2010, 20:46

[TigerSHARC]

Всё теперь стало ясно.
Объясняю суть проблемы. необходимо узнать все гармоники , которые можно полусчить с частотой дискретизации 6400. Гармоники кратны основной частоте (которая меняется от 42.5 до 57.5). Я подаю всего лишь три - такак это просто модель, что бы суть пока уловить.
Дело вот в чём - раньше я полагал, что выборка должна быть фиксирована по количеству отсчётов (жестко 128 например). Теперь понятно что происходит так :
1) Получаем частоту (fосн.), по некоторому алгоритму (давайте договоримся что частота известна с абсолютной точностью, так как частотный алгоритм - это отдельный разговор).
2) затем интерполятор на основе данных о частоте забирает выборку АЦП из кольцевого буфера (о буфере ниже). Причём забирает столько выборок сколько ему надо (в зависимости от частоты сигнала) - если частота изменяется от 42.5 до 57.5, то кол выборок может быть от 111 до 150.
3) Затем интерполятор выдаёт ровно 128 выборок (всегда!) на новой частоте дискретизации равной 128*fосн.
4) по этим выборкам получаем спектр.
_____________________________________

КОЛЬЦЕВОЙ БУФЕР

теперь главный вопрос:

Данные с АЦП поступают в кольцевой буфер. Буфер заполняется постоянно со скоростью, которая зависит от периода дискретизации (в моём случае 1/6400).
Допустим частота сигнала найдена и равна 43Гц. Интерполятору требуется 148 выборок для данного временного окна и он забирает первые 148 выборок из буфера, 149 выборка уже является первой для следующего временного окна, потом частота изменилась и скажем равна 45 Гц, тогда интерполятор берёт по 142 выборки на временное окно и так далее...
Кольцевой буфер должен быть такого размера что бы хранит предысторию сигнала за некоторое время и должен выдавать её при запросе.
Буфер заполняется постоянно и циклически. так вот проблема в том что скорее всего буфер заполняется быстрее, чем выборки извлекаются из него. Получается что, какого бы размера не был буфер, рано или поздно вершина буфера настигнет его хвост(где храняться полезные данные)...

как этого избежать при работе с кольцевым буфером.

[Andron_]

Получается что, какого бы размера не был буфер, рано или поздно вершина буфера настигнет его хвост(где храняться полезные данные)...

как этого избежать при работе с кольцевым буфером.

взять кольцевой буфер бесконечного размера, применив бесконечное число микросхем памяти.

[thermit]

TigerSHARC:
уфер заполняется постоянно и циклически. так вот проблема в том что скорее всего буфер заполняется быстрее, чем выборки извлекаются из него. Получается что, какого бы размера не был буфер, рано или поздно вершина буфера настигнет его хвост(где храняться полезные данные)...

как этого избежать при работе с кольцевым буфером.

Вопрос из области "как сделать так, чтоб 2*2 равнялось бы 5?"
Ответ очевиден - никак.

Чтобы буфер не переполнялся нужно забирать данные из него с той же средней скоростью, с которой они в буфер кладутся. Храните тада передискретизированные данные...

[TigerSHARC]

... а если делать так:
брать данные из буфера и выполнять вычисления с ними тогда когда в буфере накопились данные... дело в том что процессор успевает вычислить БПФ.
Получается так:
интерполятор запросил выборку (от 120 до 150 отсчётов), выполнилось БПФ, а пока выполнялось - копилась новая выборка, затем опять запрос выборки (от 120 до 150 отсчётов), и снова БПФ, и снова копиться выборка пока делаем БПФ...
Можно ли так сделать, если выборка берётся не фиксированной длинной (от 120 до 150 отсчётов)...?

Да это относиться не только к моему случаю... как вообще делается БПФ в реальном времени, там же тоже кольцевой буфер... и как там решена проблема с кольцевиком?

Сообщение отредактировал TigerSHARC - Feb 25 2010, 12:42

[thermit]

Если ориетироваться на число отсчетов связанное с частотой 1-й гармоники - то никаких проблем не будет.
Советую почитать что-нибудь про системы автоматической подстройки частоты. Именно она у Вас и получается.

[TigerSHARC]

Если ориетироваться на число отсчетов связанное с частотой 1-й гармоники - то никаких проблем не будет.
Советую почитать что-нибудь про системы автоматической подстройки частоты. Именно она у Вас и получается.

Если бы можно было ФАПЧ использовать аппаратно... то проблем бы не было...
А тут так альтернатива - программная подстройка...

Короче надо как то проблему с циклическим буфером решить... и непонятно причем тут 1-я гармоника (выборка есть выборка) берём от 120 до 150 отсчётов в зависимости от обстоятельств(хоть там 50 гармоник).

Может так: БПФ успевает выполниться, скажем к тому времени пришло 90 отсчётов след. окна. Тогда ждём пока не будет нужного количесвта выборок, которые "запрашивает" интерполятор. Выборка пришла - опять БПФ... итак по кругу...
не выйдет так?

Причём "окна " не должны перекрываться. Каждая выборка участвует в БПФ только со своей группой выборок (это так для справки)

Сообщение отредактировал TigerSHARC - Feb 25 2010, 12:55

[thermit]

Тогда ждём пока не будет нужного количесвта выборок, которые "запрашивает" интерполятор. Выборка пришла - опять БПФ... итак по кругу...
не выйдет так?

Выйдет.

[TigerSHARC]

Подведём итог в теме.
Предисркретизация осужествляется дробным интерполятором. Который в свою очередь "решает" какое количество выборок с АЦП взять, в зависимости от частоты основной гармоники. Далее производиться БПФ сигнала.
При этом выборки с АЦП циклически заносятся в кольцевой буфер, корректная работа которого возможна в том случае, если БПФ успеет выполнится до прихода новой выборки (при этом исключается эффект когда новая выборка может когда-то записаться в "рабочую" область).

Если всё верно - всем спасибо.

________________________________

а как ещё кроме кольцевого буфера можно решить эту проблему?

Сообщение отредактировал TigerSHARC - Feb 25 2010, 19:31

[EvgenyNik]

1. Время обработки всегда примерно одинаково, т.к. функции ресемплинга в любом случае надо сформировать 128 точек, выполнив 128*M математических операций (M - число математических операций на формирование одной выходной точки).
2. Время выполнения БПФ одинаково всегда
3. Если пропуски в наблюдении недопустимы, то суммарное время п.1 и п.2 не должно превышать время набора исходных данных для максимальной частоты основной гармоники, т.е. не должно быть больше, чем время набора 111 точек на частоте 6400
4. Если п.3 выполняется, то глубина буфера не должна быть меньше, чем число точек на минимальной частоте основной гармоники (150 точек) + число точек на максимальной частоте (111 точек), тогда пока выполняется обработка набранного окна, новые данные не затрут его.

[TigerSHARC]

Да, нужно ещё учесть что планируется хранение предыстории сигнала за некоторый период.

А вообще всё здорово получается! Спасибо всем. Думаю тема исчерпана.

[megajohn]

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

Отлично ! это то что искал - работает. Хочется добавить, если кому хочется отказатся от массивов t1 и t0, то заменить их где можно а строку double x = (t1[i]-t0[j])*Fs0-2; заменить на x = (double)i * Fs0 / Fs1 - j - 2;

[megajohn]

хоть тема и древняя, но нужная всегда.
В аттаче пример на студии10 передескритизации WAV-файла на основе исходника от bahurin`а
Когда наглядно, и можно по шагам пройтись по алгоритму, то будет проще понять.

Прикрепленные файлы

 pc_resample.rar ( 334.2 килобайт ) Кол-во скачиваний: 37