DSPSYSTEM Теория и практика цифровой обработки сигналов

Информация о пользователе

Привет, Гость! Войдите или зарегистрируйтесь.



фильтр фарроу

Сообщений 1 страница 20 из 20

1

Хочу передискретизировать сигнал.
Чистый синус частотой 45.35 гц. Исходная частота дискретизации 6400 Гц. Выборок 128. Хочу передискретизировать так, чтобы во временном окне помещалсяровно период сигнала. Как это сделать?
Использую код, представленный на вашем сайте:

#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 = 6400;            // исходная частота дискретизации
    double Fs1 = 5804.8;            // частота дискретизации после ресэмплинга

    int N = 128;// количество исходных отсчетов                   
    int K = (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] = sin(2*M_PI*45.35*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;
}

Что то не очень получается(((
Где я ошибаюсь?

P.S. Спасибо за сайт!

0

2

один период синусоиды с частотой 45.35 Гц составляет 1/45.35 = 0,02205 сек. При этом 128 отсчетов с частотой дискретизации 6400 составляют во времени 128/6400 = 0,02 сек, а 128 отсчетов с частотой дискретизации 5804.8 Гц составляет во времени: 128/5804.8 = 0,02205. Таким образом, для того чтобы получить полный период на выходе ресэмплера должно быть 128 отсчетов. Поскольку частота понижается, то необходимо чтобы на входе  было N = 128*6400/5804.8 = 142 отсчета. Тогда для вашего случая:

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

0

3

Справедливо ли такое преобразование? По сути получаем просто вырезанный период сигнала. кстати - хороший метод борьбы с растеканием спектра. Если конечно сигнал периодичен.

0

4

а что собственно смущает?

0

5

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

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

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

0

6

чтото я не понял это что за монолог был?  %-)

0

7

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

так вот, к сожалению буду голословен, но должно быть примерно так:
"...интерполятор работает так что берет такое количество отсчетов что в них содержится ровно один период, независимо от частоты "(это цитата специалиста)
И это действительно так. Реализуется с помощью фарроу-фильтра. (это уже моё)

незнаю как пока что... (тоже моё)

0

8

если частота известна то пример выше показывает как это сделать.

0

9

Да, но как же быть с тем что в таком случае часть сигнала фактически просто отбрасывается. или я не прав?
Интерполятор должен не вырезать один период, а брать один период. Это означает, что все последующие выборки идут в буфер следующего БПФ.
А у здесь просто отбрасывается.

P.S. Да, кстати, спасибо за сайт. Много бесценных вещей!

Отредактировано Franky (2010-02-24 16:27:26)

0

10

ну почему же отбрасывается? Если вы взяли первые 142  это на частоте 6400 во времени составляет 0,0221875 сек. Из них в результате ресэмплинга получили 128 отсчетов с частотой  5804.8 Гц. Во времени это 0,02205 секунды. Таким образом интервал обработки отличается на 0.1375 мс. После следующие 142 отсечета снова превратятся в 128 и т.д. вы ничего не отбрасываете просто пересчитываете время. Естественно, что между интервалами обработки будет дырка в 0.1375 мс, что соответствует фазовому шуму 0.3 градуса при частоте сигнала 45 Гц.

0

11

например сигнал имеет частоту 55 Гц.
Мы берём фиксировано 150 отсчётов (это априори) на частоте 6400 (чтобы поймать период самой низкой из допустимых частот синуса 42,5Гц). Тогда размер временного окна 0,0235 с.
Синус частотой 55Гц содержит период 0,0182с.
делаем передискретизацию на частоту (128*55) = 7040 Гц и 128 отсчётов. Получаем строго период.
Вот код матлаб:

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));

При этом весь сигнал, что в интервале 0,01818...0,0235с. просто отбрасывается... а разрывов между окнами анализа быть не должно...

Неужели нельзя просто взять исходные 128 отсчётов на фиксированном временном окне и просто подогнать под них сигнал...(((

Отредактировано Franky (2010-02-24 17:34:52)

0

12

Просто вы никак не уловите суть. Окно должно меняться в зависимости от частоты. Для низкой частоты вы должны входные данные обрабатывать кусами по 150 отсчетов, а для самой высокой по 120 например. И так для кажой частоты из сетки. Тогда вы всегда будете обрабатывать стык в стык ничего не теряя.

0

13

С вашей помощью постараусь уловить суть.
Скажите, а разве можно брать разное количесвто отсчётов ??? Я всегда думал, что АЦП выдаёт определённое количесвто отсчётов в заданном временном окне с определённой частотой дискретизации...
P.S. было бу удобнее пообщаться в ICQ.
Спасибо за помощь.

Отредактировано Franky (2010-02-24 18:30:10)

0

14

да но можно ведь взять 128 отсчетов с этого ацп и временное окно будет одно,а если взять 150 то временное окно будет другое. Потом можно из 150 сделать 128 и получите что отсчетов 128 (как и в первом случае) но окно временное другое, так как изначально было 150 отсчетов с ацп

0

15

дело в том ч то если брать 150 отсчётов всегда, то получиться что часть сигнала вырезана и отброшена (в моём случае максимальное максимальный интервал этого отбрасываемого куска 6 мкс. поскольку диапазон частот синуса 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 выборки на временное окно и так далее...
Кольцевой буфер должен быть такого размера что бы хранит  предысторию сигнала за некоторое время и должен выдавать её при запросе.
Буфер заполняется постоянно и циклически. так вот проблема в том что скорее всего буфер заполняется быстрее, чем выборки извлекаются из него. Получается что, какого бы размера не был буфер, рано или поздно вершина буфера настигнет его хвост(где храняться полезные данные)...

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

Отредактировано Franky (2010-02-25 15:11:50)

0

16

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

0

17

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

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

0

18

Franky написал(а):

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

да тогда все хорошо

0

19

Чтобы не плодить темы задам вопрос по Фарроу тут.

Интересно, а в чём разнится когда мы передискретизируем фильтром Фарроу с 25 кГц на 5кГц, и когда мы передискретизируем с 250кГц на 5 кГц. Иначе говоря, чем будут отличаться два сигнла на выходе фильтра фарроу (оба с частотами 5 кГц) в зависимости от степени передисретизации?

0

20

ничем в данном случае просто децимация

0