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

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

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

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

тогда, как я понимаю, 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

Здравствуйте!
Может быть вот такой вариант Вас устроит:
% 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

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

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

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

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

За код - пожалуйста, тем более это ваш код в общем то с небольшими изменениями.
По полифазным КИХ я выше написал неверно - вполне их тоже в этой ситуации можно использовать.
Извините, 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

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

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

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

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

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

дело в том что нужно обеспечить отсутсвтие растекания спектра, если основная гармоника колеблется в пределах 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 - Feb 24 2010, 12:33

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

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

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

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

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

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

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

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