fft/ifft для последовательности длиной 768 (=512+256) Опции

[PriBoris]

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

Сообщение отредактировал PriBoris - Jul 27 2010, 08:49

[bahurin]

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

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

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

данные поступают в реальном времени блоками по 512.

Сообщение отредактировал bahurin - Jul 27 2010, 10:07

[Xenia]

Кажется Вы пытаетесь сказать новое слово в цифровой обработке сигналов...

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

[bahurin]

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

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

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

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

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

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

Сообщение отредактировал bahurin - Jul 27 2010, 10:47

[Xenia]

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

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

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

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

[bahurin]

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

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

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

Код

Fs = 100;       %частота дискретизации
f0 = 20.1;      %частота сигнала
N = 512;        %количество отсчетов
t = (0:N-1)/Fs; % время

s = sin(2*pi*f0*t); % исходный сигнал 512 отсчетов
S = 20*log10(abs(fft(s)));  % спектр исходного сигнала без всякого сглаживания как есть
s0 = linspace(s(N),s(1),N); % заполняю линейный тренд от последнего отсчета до первого
s1 = [s,s0];                % добавляю к исходному сигналу линейный трэнд
S1 = 20*log10(abs(fft(s1)));    %спектр сигнала с добавлением 1024 отсчета
freq = (0:N-1)*Fs/N;            %частота при 512 отсчетах
freq1 = (0:2*N-1)*Fs/(2*N);     % частота при 1024 отсчетах

subplot(211),plot(1:2*N, s1, 1:N,s,'r'), grid on;      % вывожу на график
subplot(212),plot(freq1,S1,freq,S,'r'), grid on;      % вывожу на график

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



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