Какой эффективный алгоритм для длин такого типа (2^n+2^(n-1)) выбрать ?
Вопрос возник при попытке реализации фильтрации с помощью быстрой линейной свёртки. Длина фильтра почти 256 точек, данные поступают в реальном времени блоками по 512.

Быть может здесь было бы проще дополнить массив 768 нулями до массива 1024 и действовать обычным образом?

Не совсем понятно что Вы делаете?
Насколько я смог понять (а если я понял не правильно, то поправьте меня), Вы пытаетесь отфильтровать сигнал переведя его в частотную область при помощи БПФ и после фильтрации переведя сигнал из частотной области обратно во временную при помощи обратного БПФ. Для этого Вам необходимо сделать БПФ длиной 768 точек.

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

Да,я так и хочу сделать. Просто на всякий случай спросил, вдруг есть возможность сэкономить.

...это ерунда, причем полная. Вон товарищ связист дело говорит...

768 = 3*2^8
последняя стадия по основанию 3, остальные - по основанию 2. или 16.
C перестановками, правда будет возня. Может, действительно проще будет взять на 1к точек?

Если вы достигните своего желания получить FFT на 768-точечном массиве, то у вас получится циклическая свертка. А вам оно надо? В большинстве практических случаев нужна простая свертка, а не циклическая, поскольку образ, находящийся в массиве 768, вряд ли периодический.
Дополнение же нулями до 1024 здесь очень удачно еще и тем, что как раз добавляет к исходному массиву данных нулевой кусок той же самой длины, как и длина функции, с которой станут сворачивать. Такое добавление является как раз минимальным для того, чтобы вместо циклической свертки получить нормальную.
Более того, даже если бы ваш массив был длиной 512, что позволило бы легко преобразовать его в FFT-образ, то и в том случае стоило бы подумать о целесообразности дополнения его нулями до 1024, чтобы избавиться от циклической свертки. А у вас длина массива просто идеальна для получения нормальной свертки, т.к. 768+256=1024.

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

Цикличность свертки это и есть причина краевых эффектов

если надо поставить фильтр, например ФНЧ, то быстрее всего это БИХ - фильтр рассчитать и использовать. В этом случае выши 256 точек сведутся к БИХ фильтру порядка не больше 10.

Если же очень хочется ких фильтр, то тогда 512 точек сигнала переводите в спектр и потом уже спектр умножаете на частотную характеристику фильтра, и переводите через ifft. Как-то я не нашел в вашем вопросе 738 точек, поскольку:

Попытка не пытка! Тогда уж я еще одно новое слово скажу
В чем, собственно, недостаток FFT в практическом смысле? В прежде всего в его циклическом характере. FFT преобразует массив так, как будто он закольцован, хотя на самом деле этого нет.
Вот практический тому пример. Положим, мы преобразуем массив данных, полученных от АЦП. Причем измеренный сигнал имеет некий линейный тренд, выражающийся в том, что уровень сигнала в первой точке массива не совпадает с уровнем в конечной (например, база сигнала постоянно повышается со временем). Такое поведение входного сигнала довольно обычно и называется дрейфом.
Однако FFT-алгоритм нам на зло станет рассматривать этот массив, как свернутый в кольцо, а потому усмотрит в нем чудовичный разрыв непрерывности в том месте, где происходит склейка в кольцо. Из этого "уступа" в частотном образе сигнала появятся большие вклады не только высокочастностных составляющих, но и низкочастотных. И все лишь для того, чтобы этот уступ воспроизвести. Можно, конечно, окнами с ними бороться - это традиционный путь. А вот я "изобрела" совсем простой путь, но исключительно эффективнный. Он состоит в том, что надо дополнить массив (примерно в два раза до следующего размера, кратного степени числа 2), но только не нулями! Вместо нулей следует в этом месте провести НАКЛОННУЮ ПРЯМУЮ ЛИНИЮ, которая одним свои концом исходит из последней точки исходных данных, а другим концом должна дойти до конца рассширенного массива именно на тот уровнь, соотвествующий первой точке исходных данных. Тогда при склейке в кольцо окажется, что переход получился максимально плавным.
Такой метод практически не дает паразитных вкладов в частотный спектр, возникающих из-за перепада уровней на концах массива данных, а потому и не нуждается ни в каких окнах или дополнительных фильтрах.

Для вычисления фильтрованного значения первой точки свежего блока используются 255 (256-1) точек предыдущего блока.

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

2. Если уж вы хотите сказать свое слово, то это слово должно быть обоснованным. А так получается что вы усилием воли решили а давайте сделаем линейную интерполяцию. А почему собственно линейную? А может лучше будет если кубическим полиномом соединить, тогда можно еще и производные выровнять и устранить точки перелома? А если полином 5 степени или 7-ой? И почему вы удваиваете длину выборки, а не увеличиваете ее в 4 раза например или в 8? Короче не убедительно это все. Кроме того, я бы еще проанализировал, что будет если просто удвоить количество отсчетов и наложить сверху весовое окно? Мне кажется спектр чище будет, чем с вашей интерполяцией.

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


если у вас предыдущие 512 отсчетов и следующие 512 во времени не разрываются, тогда бих фильтр вам предпочтительнее. Если они имеют дырку между собой, то брать предыдущие отсчеты нельзя. Если же все таки хотите через fft, то вам необходимо применить обработку с перекрытием. тогда надо делать fft на 512 отсчетов, умножение на чх фильтра, затем ifft, потом сдвиг на 256 отсчетов и снова fft на 512. Т.о. вы получите аналог ких фильтра, но суть в том, что надо делать fft на удвоенную длительность всегда, чтобы было перекрытие.

Т.е. проблему решит не какое-то дополнение нулями - ну же, Xenia, у меня блок где замерян постоянный ток. 11111...111 (V) - как его не дополняй - спектр станет заметно отличным от простой палки на 0-ой частоте... Или вон пример который привел товарищ bahurin - тоже очевиден для мыслеэксперимента.
Проблему решит:
1) Выбор других базисных функций и метода построения базиса.
2) Выбор "иного" математического мира - и так между прочим тоже делают...

Остальное прочел... Здорово. Но это частные случаи.

P.S.: был задан конкретный вопрос и товарищ связист дал на него самый компетентный ответ - смысле на самый первый вопрос(уточнение аФФтАр топика дал просто чумовое ). Дальше понеслось - и дополнение 0-ми, и IIR... Сейчас пойдут теоретико-числовые алгоритмы и кольца Галуа...

P.P.S.: Xenia, перепишите DCT(FCT) для своего mp3 плеера или JPEG просмотрщика - так "грамотно". С "не кольцевой" сверткой. Послушайте(посмотрите) его минут 10. А выводы опубликуйте на форуме...

P.P.P.S.:
Мда... Это диагноз... PriBoris, без обид, но Вы жжоте...

Почти это как? 255.5?

делайте по 2-а FFT по 256 точек...
или одно по 512 на блок...

Физический смысл быстрой свертки в том, что у Вашего фильтра есть определенные АЧХ/ФЧХ(которое однозначно связано с его "импульсной"). Мы определяем его и накладываем операцией умножения(свертка-фильтрация над временным рядом соответствует умножению в частотном отображении) на частотный спектр сигнала. Размер блока может быть любым - впринципе, даже меньше длинны фильтра - но тогда его нужно будет пересемплировать с другой частотой Ничего из уже обработанного блока хранить не надо. Это не FIR. А все ускорение за счет того, что переход в частотную область с помощью FFT, посемпловое умножение, и обратный переход уже где-то начиная с 128 точек требует меньше операций чем свертка "в лоб".

Кстати блоки по 256 точек, 32-х битная система и тип данных int(если у Вас, PriBoris, так) - это и в самом деле лучше на кольцах Галуа. Зачем нам комплексность и округленные значения sin(t)/cos(t) !?

В чем вопрос?

Где ж это в практических случаях бывает, чтобы изменяемая гармоника да уложилась бы что в целое число периодов в массив данных? Вы привели выгодный для себя пример . А давайте я вам приведу пример, когда гармоника (предположим, косинусоида) заканчивается "неудачно". Скажем начала косинусоида значение +1, а занончилась -1 (минус один), пол периода ей не хватило. Знаете, как гадость тогда получается? Одной палкой даже не пахнет . А вот мой метод с такими случаями отлично справляется. И пусть там не чистая палка, но боковых липествов посчи нет.
И что из того, что нелинейно? Откуда у вас такая приверженность к линейности и отвращение к нелинейности? Линейным преобразованием из сигнала даже большую гадость можно сотворить, чем линейным . Опять же сейчас речь идет не об усилительных каскадах, где нелинейные искажения почитаются грехом, а о совершенно иной задаче.


Основания моих слов имеются, причем они примерно такие же, как причины выбора формы окна при фильтровании в частотной области. Ведь никто же не режет отсечкой с некоторой частоты, поскольку знает, что за этим поледует. А последует "зашумленность" той самой граничной частой (при обратном FFT), по которй обрезали. Поэтому если и давят частоты, то осторожно - так, чтобы коэффициент подавления сходил к нулю плавно. И тут кто во что горазд: и треугольное окно, и косинусоидное, и гауссовый профиль, и т.п.
Дополнительный треугольник, достоенный в исходной области данных - наименьшее зло, по сравнению со всеми остальными случаями, когда краевые точки массива сильно отличаются по величине. Причем сам треугольник ведет себя на редкость смирно. Можно построить "домик" (линейное возрастание до середины массива с убыванием его во второй половине массива к начальному значению). Тут "мусор", порождаемый таким треугольником слаб, а амплитуды частот изменяются непрерывно, без каких либо максимумов. Анализировать спектральные данные (наложенные на этот треугольник) это не мешает. Опять же естественные процессы линейного дрейфа сигнала есть ни что иное, как добавление аналогичного трегольника. Поэтому такую форму искажения тоже можно считать естественной, присущей реальным задачам обсчета временных трендов.

По большому счету - нелинейность необратима. Помню Tanya каГ-бЕ намекала нам, когда я сам интересовался этим вопросом. Вернее обратима если точно знать характеристику нелинейного звена и первоначальный входной сигнал. Но Хевисайд уже смотрит на Вас с укоризной...

все это эмпирика, рассуждения. Не первый раз сталкиваюсь с ситуацией когда человек что-то якобы "изобрел" и пытается обосновать. То что вы предлагаете устранить скачок фазы при вычислении fft это понятно. Подобным образом умножение на окно также приводит начало сигнала к его концу. Вы делаете данное приведение интерполяцией полиномом первой степени. Так вот прежде чем говорить о своем достижении, неплохо было бы провести сравнительный анализ своего алгоритма и существующих. Например насколько ваш метод давит боковые лепестки при растекании спектра? Я почему-то думаю, что это будет зависеть от формы сигнала. Для какого-то сигнала будет одно подавление, а для другого - другое. Вот именно в этом и есть нелинейность, когда подавая различные исходные данные вы получаете различный результат. Вот например если у меня так получилось, что частота дискретизации почти кратна частоте синусоиды, то я бы и без вашего метода имел бы отличный спектр, а применив ваш метод я фигню полную получу, потому что вместо синусоиды у меня будет синусоидальный импульс и как следствие - уровень боковых - 13 дБ. В этом и есть нелинейность вашего метода.

а вот вам код матлабовский для построения спектра:


А вот вам спектр синий - ваш метод, красный вообще без ничего как есть

Нажмите для просмотра прикрепленного файла

синий это ваш метод. Мало того что лепестки выросли, еще и постоянка появилась из-за вашего изобретательства.

Спасибо,с алгоритмом Кули-Тьюки разобрался. В матлабе прокрутил вариант 3*256, всё сходится.


Нет,имелось ввиду,что длина фильтра после различных доработок надфилем получается в районе 210...230.
256 - это я взял с запасом на будущее,округлив до ближайшего красивого (для хранения и использования) числа



Обясните пожалуйста, я не понимаю что я по Вашему мнению планирую делать не так.
Мне нужно реализовать FIR длиной 256, который я собственно реализовал в лоб. Получается затратно по ресурсам.
Я планирую сделать так:
Пришел новый набор данных из 512 точек.
Я создаю буффер их 768 точек, в него в позиции с 1 по 256 я копирую последние 256 точек из предыдущего набора данных, а в позиции с 257 по 768 я копирую свежие 512.
Делаю fft длиной 768.
Домножаю результат на АЧХ/ФЧХ моего фильтра, дополненного нулями до длины 768.
Полученное произведение пропускаю через ifft.
Точки с с 1 по 256 выкидываю, а точки с 257 по 768 считаю результатом фильтрации свежего блока данных.

Лучше делать два таких прохода на блок, но по 512 точек ?
Или один проход, но избыточный, на 1024 точек, чтобы не заморачиваться с Кули-Тьюки (в котором я возможно потреряю много времени на накладных расходах при переиндексации и т.д.)?

Ваше требование некорректно . В вашем примере FFT применяется для сжатия сигнала (к более компактной форме представления/хранения), а потому о концевых эффектах там просто нет речи - важно лишь, чтобы обратное FFT воспроизводило исходный сигнал достаточно хорошо.
Здесь мы имеем тот эффект, что последовательное применение FFT и iFFT (если ничего не обрезать в промежуточном частном образе) тождественно воспроизведет (дискретный) оригинал. И это несмотря на любые краевые эффекты.
Но если нам нужно частотное представление само по себе, например, для гармонического анализа анализа или получения свертки, то с "краевыми эффектами" мириться нельзя. Ибо разрыв (ступенька между уровнями) на краях массива порождает при FFT дополнительные частоты, которые на самом деле во входном сигнале отсутствуют. И в этом можно легко убедиться, если использовать обычный метод FT (не Fast).
В принципе можно с определенностью сказать следующее: все три метода (зацикливание, дополнение нулями и дополнение линейным трендом) грешат против истины. В самом деле, при вычислении сверки прямым способом в конце массива создается ситуация, когда в результате относительного сдвига двух сворачиваемых функций приводит к неопределенной ситуации - первая функция "кончилась", выйдя за область, где мы знаем ее значения, а вторая функция еще нет. Если в этой ситуации мы решаем довольствоваться частичной суммой произведений, то это будет эквивалент заполнения недостающей части нулями. А если же решим, что "закончившаяся" функция периодична и решаем скопировать недостающие точки из ее начала - это то, что делает FFT. А мой "трегольник" - это попытка обмануть FFT, подсунув ей данные, на которых принудительная "циклизация", которую делает FFT, приносит наименьший вред (в случаях, когда периодичность исходной функции заведомо отсутствует).
Величина этого вреда может колебаться в широких пределах, в зависимости от того, насколько совпадают все эти три предположения (о дальнейшних значениях функции) с реальностью. В случае, если исходный сигнал - чистая гармоника, да еще и уложившаяся в интервал целым числом периодов, то, конечно же, FFT тут победит, т.к. ее "преположение" о перирдическом характере сигнала полностью оправдалось. Да вот только не бывает на практике таких удачных совпадений. И если вдруг окажется, что за переделами участка сигнал какой-то другой, то преположение о периодичности и цельно-кратности периодов окажется лажей. В этом случае ошибка может быть очень большой. Заполнение нулями - это согласие на среднюю ошибку. Она окажется больше, если в периолд удалось вписаться, или меньше, если вписаться не удалось. Т.е. здесь мы уповаем на нулевое среднее того сигнала, продолжения которого не знаем. Вообще-то это не такое уж плохое предположение.
Ситуация становится из рук вон плохой, когда мы решились на нулевое среднее, но продолжаем использовать метод FFT. Вот тут-то нас и подстерпегают сюрпризы из-за того что FFT-пребразование апроксимирует гармониками те разрывы непрерывности, которые имеют место (в точке перехода крайней точки на нуль и точку перехода нуля на первую точку). И вот именно с этим эффектом борется "метод треугольника", делая этот переход наиболее плавным.

А почему не 512?


А зачем, простите? Это Вы так перекрытие делаете или что?


???


Тут зависит от задачи - иногда нужно дополнить, а иногда масштабировать...


?????????????????????


Это мегаерундище...

Еще раз - это не FIR. Совсем не FIR. Т.е. о свертке забываем. А значит ничего не храним. Теперь у нас только умножение в частотной области - блочный алгоритм.




Лучше сесть и спокойно разобраться - как же оно - правильно.

Вы все верно спланировали. Это метод перекрытия с накоплением.

Я бы с удовольствием, но реализовав в матлабе описанное, вижу что результат этого метода не отличается от результата встроенной функции filter(b,1,in), которой я скармливаю свой оригинальный фильтр длиной 219.
В иностранной литературе метод называют Overlap-Save. (см.вложение)

Золотые слова, юрийвенедиктович...

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


Йа кажеЦЦо понял причину Ваших душевных терзаний . Так вот - сядьте и немного повоображайте. И подумайте когда произойдет выделенное мной если один из сигналов бесконечен. А в случае RT данных это очень удачная идеализация.


Да, но это будет в конце рабочего дня когда оператор выключит приборы.


...ну, мягко говоря, FFT не так работает.


Обманывать - непедагогично!


Это черная магия от ЦОС.


...ну раз в иностранной - тода, да... Я умываю руки, собственно кто йа такой чтобы идти против лавины иностранной науки!!!

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

Никаких 512-и быть тут не может, т к результат свертки 2-х последовательностей длин 512 и 256 будет иметь длину 767. Т е размер дпф должен иметь длину минимум 767. Ближайшая внятная размерность 768. Что касается 1024-х точек - уверен, что дпф на 768 точек будет быстрее. Это справедливо, если нет аппаратной поддержки перестановок. (Например, если есть аппаратная поддержка битреверсной адресации, которая используется в алгоритме к-т по основанию 2, то есть выигрыш в числе операций). Так что не парьтесь, Вы на правильном пути.

ps
Можно манипулировать дпф сигнала используя окна без умножения на дпф фильтра, как пытались посоветовать некоторые товарищи. Но это будет совсем другая история, и размер 512 тут далеко не оптимален... Хотя, все зависит от требований к точности.

Под двумя проходами по 512 я имел ввиду :
первый проход : 512 - состоит из второй половины предыдущего блока и из первой половины нового блока данных
второй проход : 512 - состоит из первой и второй половин нового блока данных
то есть на каждый блок данных длиной 512 я прохожусь 2 раза сверткой с длиной 512, и после каждой свертки получаю 256 фильтрованных значений

Сегодня этот вариант реализовал и проверил в процессоре. Выигрыш по сравнению с прамой реализацией FIR составил 6 раз. На данном этапе меня это устраивает. Если в будущем не будет хватать, то буду делать 768-точечное FFT.
Большое спасибо всем откликнувшимся!

Согласно самому определению свертки двух конечных рядов во временном представлении - да. Но товарищу нужно наложить фильтр. Просто фильтр наложить!!! И хрен он его наложит на бесконечный ряд согласно этому определению. Классическая длинна блока, классический алгоритм БПФ, делается в течении дня вместе с тестами... Зачем вся эта ерунда!!!??? Разрабатывать БПФ странной длинны... Добавляется Х.з. сколько операций. В крайнем случае фильтр можно немного "предисказить" для улучшения результата. Хотя - мне за это не платят, так что чего это я ...

А доказательство своей правоты наличием невнятной англоязычной статьи - это супер.

КАроче - школота и дилетанты.

Что он добросовестно и проделал.
КИХ фильтрация и апериодическая дискретная свертка - одно и тоже.
Изобретать тут нечего. Все уже изобретено до нас.

зы:
На счет школоты и дилетантов - эт Вы погорячились...

Да. Но надо различить 2 случая. 1. имеется массив данных s длины N. Выход ких фильтра с заданной импульсной характеристикой h равен линейной свертке s*h. Но есть еще один случай, когда длина массива s неизвестна. Например данные идут кусками по 512 отсчетов. В этом необходимо чтобы при поступлении следующего куска данных, ячейки памяти ких фильтра хранили информацию о предыдущем куске. И если вы просто в лоб выполните линейную свертку каждого куска и потом их состыкуете, то результат будет отличаться от того, что вы получите если выполните линейную свертку всех кусков сразу. Для того чтобы хранить данные о предыдущем куске делают обработку с перекрытием (overlapped), т.е. каждый раз при поступлении следующего куска захватывают часть предыдущего чтобы заполнить ячейки памяти фильтра.

А какого размера относительно длинны массива должно быть перекрытие?

Длине импульсной характеристики без 1.



Это круговая (циклическая) свертка.



А это вычисление апериодической свертки через круговую. Вопщем все верно...