Справочник функций

Ваш аккаунт

Войти через: 
Забыли пароль?
Регистрация
Информацию о новых материалах можно получать и без регистрации:

Почтовая рассылка

Подписчиков: -1
Последний выпуск: 19.06.2015

Фильтрация методом преобразования Фурье

505
17 февраля 2007 года
vAC
343 / / 28.02.2006
Возникла проблемка:
В программе есть фильтр, схема работы следующая образом:
1. БПФ сигнала
2. memcpy(&pSignal[FilterAperture/2], pSignal[Aperture-FilterAperture/2], FilterAperture/2*sizeof(COMPLEX_FLOAT32));
Здесь происходит "укорачивание" спектра
4. Взвешивание "укороченного" спектра окном Кайзера
3. ОБПФ "укороченного" спектра

Aperture - длина сигнала
FilterAperture - длина "укороченного" сигнала

А проблема в том, что получается сигнал, в котором первый отсчет "темнее" остальных. Т.е. при пропускании изображения строками, получается темная полоса вначале...
Вся фигня в пункте 2 и соответственно 3, т.е. даже без взвешивания получается затемнение. Просто теоретического обоснования такого "укорачивания" я не нашел.
Точнее проблему сформулировать так:
как из частотной характеристики из N отсчетов получить характеристику состоящей M отсчетов, M<N, чтобы при обратном преобразовании небыло никаких затемнений...
551
17 февраля 2007 года
Pavia
357 / / 22.04.2004
2. Незнаю зачем тебе понадобилось выкидовать часть спектора.
memcpy(&pSignal[FilterAperture/2], &pSignal[Aperture-FilterAperture/2], FilterAperture/2*sizeof(COMPLEX_FLOAT32));
&-пропущен

4. Взвешивание "укороченного" спектра окном Кайзера
3. ОБПФ "укороченного" спектра

Порядок цифр и пунктов не совподает.
Чето не понятно как-то не к месту применяется окно Кайзера.

Насчет затемнения, ты когда выкидываешь часть спектора, с ним уходит энергия. Но не в этом суть, затемнение тогда было бы на всех отчетах.
505
17 февраля 2007 года
vAC
343 / / 28.02.2006
Цитата: Pavia
2. Незнаю зачем тебе понадобилось выкидовать часть спектора.
memcpy(&pSignal[FilterAperture/2], &pSignal[Aperture-FilterAperture/2], FilterAperture/2*sizeof(COMPLEX_FLOAT32));
&-пропущен

4. Взвешивание "укороченного" спектра окном Кайзера
3. ОБПФ "укороченного" спектра

Порядок цифр и пунктов не совподает.
Чето не понятно как-то не к месту применяется окно Кайзера.

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



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

bool Filter_HFF::TransformCol(int iCol, int Height, class FilterBuffer *pIn, class FilterBuffer *pOut)
{
int i;

pIn->CopyCol(m_Signal, iCol, 0, Height);
memset(((COMPLEX_FLOAT32*)m_Signal)+Height, 0, (m_Aperture-Height)*sizeof(COMPLEX_FLOAT32));

m_SignalFT.FFT(m_Signal);

memcpy(&((COMPLEX_FLOAT32*)m_Signal)[m_FilterAperture/2], &((COMPLEX_FLOAT32*)m_Signal)[m_Aperture - m_FilterAperture/2], m_FilterAperture/2*sizeof(COMPLEX_FLOAT32));
memset(&((COMPLEX_FLOAT32*)m_Signal)[m_Band/2], 0, (m_FilterAperture - m_Band)*sizeof(COMPLEX_FLOAT32));

for(i = 0; i < m_Band/2; i++)
{
((COMPLEX_FLOAT32 *)m_Signal).real *= float(m_pWeight[i + m_Band/2]);
((COMPLEX_FLOAT32 *)m_Signal).imag *= float(m_pWeight[i + m_Band/2]);

((COMPLEX_FLOAT32 *)m_Signal)[m_FilterAperture - m_Band/2 + i].real *= float(m_pWeight);
((COMPLEX_FLOAT32 *)m_Signal)[m_FilterAperture - m_Band/2 + i].imag *= float(m_pWeight);
}
for (i = 0; i < m_FilterAperture; i++)
{
((COMPLEX_FLOAT32 *)m_Signal) = ((COMPLEX_FLOAT32 *)m_Signal)[i*m_Aperture/m_FilterAperture];
}
m_FilterFT.IFFT(m_Signal);

COMPLEX_FLOAT32 *pInput = m_Signal;
BYTE *pOutput = *pOut;
int Pitch = pOut->Pitch();

// pInput += m_RefLen;
pOutput += iCol*sizeof(COMPLEX_FLOAT32);
for (i=0; i<m_FilterAperture; i++)
{
((COMPLEX_FLOAT32*)pOutput)->real = sqrt(pInput->real*pInput->real+pInput->imag*pInput->imag);
((COMPLEX_FLOAT32*)pOutput)->imag = 0;
pInput++;
pOutput += Pitch;
}
return true;
}

А часть спектра выкинуто для автоматического получения укороченного сигнала при обратном Фурье (скажем в 16 раз меньшего, чем исходный). А насчет окна Кайзера это я сказать не могу зачем...так было сделано до меня, а моя задача состоит совсем в другом, сделать которую как раз и мешают эти полосы.

505
17 февраля 2007 года
vAC
343 / / 28.02.2006
А насчет энергии...как раз энергию я не выбрасываю, а наоборот, для Фурье меньшего размера такой спектр получиться с большей энергией.
Вообщим при обработке тестовой последовательности следующим образом:
CFastFT_RI4_RI4 ft,ft2;
COMPLEX_FLOAT32 *data;
int k;

data = _new COMPLEX_FLOAT32[16];
for (k=0; k<16; k++)
{
data[k].real = k;
data[k].imag = 0;
}
ft.SetFourieOrder(4);
ft2.SetFourieOrder(3);
ft.FFT((float*)data);
memmove(data+4, data+12, 4*sizeof(COMPLEX_FLOAT32));
for (k=0; k<8; k++)
{
data[k].real /= 2;
data[k].imag /= 2;
}
ft2.IFFT((float*)data);
for (k=0; k<8; k++)
{
data[k].real = sqrt(data[k].real*data[k].real+data[k].imag*data[k].imag);
data[k].imag = 0;
}
Получаю следующее преобразование:
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15
->
4.03, 1.09, 4.49, 5.82, 8.01, 10.21, 11.54, 15.03
Здесь уже на первом "светлый", на втором "темный", а дальше более или менее.
551
17 февраля 2007 года
Pavia
357 / / 22.04.2004
Цитата: vAC
Извиняюсь за опечатки, просто "упростил" исходник, а так он выглядит на самом деле:

bool Filter_HFF::TransformCol(int iCol, int Height, class FilterBuffer *pIn, class FilterBuffer *pOut)
{
int i;

pIn->CopyCol(m_Signal, iCol, 0, Height);
memset(((COMPLEX_FLOAT32*)m_Signal)+Height, 0, (m_Aperture-Height)*sizeof(COMPLEX_FLOAT32));

m_SignalFT.FFT(m_Signal);

memcpy(&((COMPLEX_FLOAT32*)m_Signal)[m_FilterAperture/2], &((COMPLEX_FLOAT32*)m_Signal)[m_Aperture - m_FilterAperture/2], m_FilterAperture/2*sizeof(COMPLEX_FLOAT32));
memset(&((COMPLEX_FLOAT32*)m_Signal)[m_Band/2], 0, (m_FilterAperture - m_Band)*sizeof(COMPLEX_FLOAT32));

for(i = 0; i < m_Band/2; i++)
{
((COMPLEX_FLOAT32 *)m_Signal).real *= float(m_pWeight[i + m_Band/2]);
((COMPLEX_FLOAT32 *)m_Signal).imag *= float(m_pWeight[i + m_Band/2]);

((COMPLEX_FLOAT32 *)m_Signal)[m_FilterAperture - m_Band/2 + i].real *= float(m_pWeight);
((COMPLEX_FLOAT32 *)m_Signal)[m_FilterAperture - m_Band/2 + i].imag *= float(m_pWeight);
}
for (i = 0; i < m_FilterAperture; i++)
{
((COMPLEX_FLOAT32 *)m_Signal) = ((COMPLEX_FLOAT32 *)m_Signal)[i*m_Aperture/m_FilterAperture];
}
m_FilterFT.IFFT(m_Signal);

COMPLEX_FLOAT32 *pInput = m_Signal;
BYTE *pOutput = *pOut;
int Pitch = pOut->Pitch();

// pInput += m_RefLen;
pOutput += iCol*sizeof(COMPLEX_FLOAT32);
for (i=0; i<m_FilterAperture; i++)
{
((COMPLEX_FLOAT32*)pOutput)->real = sqrt(pInput->real*pInput->real+pInput->imag*pInput->imag);
((COMPLEX_FLOAT32*)pOutput)->imag = 0;
pInput++;
pOutput += Pitch;
}
return true;
}

А часть спектра выкинуто для автоматического получения укороченного сигнала при обратном Фурье (скажем в 10 раз меньшего, чем исходный). А насчет окна Кайзера это я сказать не могу зачем...так было сделано до меня, а моя задача состоит совсем в другом, сделать которую как раз и мешают эти полосы.



Вроде бы тут так
pIn->CopyCol(m_Signal, iCol, 0, Height);
memset(&((COMPLEX_FLOAT32*)m_Signal)[Height], 0, (m_Aperture-Height)*sizeof(COMPLEX_FLOAT32));

А вот это похоже на то что реализован фильтр. Только как то очень каряво.
memset(&((COMPLEX_FLOAT32*)m_Signal)[m_Band/2], 0, (m_FilterAperture - m_Band)*sizeof(COMPLEX_FLOAT32)); // Почему от m_Band/2 ??? помойму надо от m_Band.

for(i = 0; i < m_Band/2; i++)
{
((COMPLEX_FLOAT32 *)m_Signal).real *= float(m_pWeight[i + m_Band/2]); // а сдесьбы я -i сделалбы
((COMPLEX_FLOAT32 *)m_Signal).imag *= float(m_pWeight[i + m_Band/2]); // а сдесьбы я -i сделалбы

((COMPLEX_FLOAT32 *)m_Signal)[m_FilterAperture - m_Band/2 + i].real *= float(m_pWeight);
((COMPLEX_FLOAT32 *)m_Signal)[m_FilterAperture - m_Band/2 + i].imag *= float(m_pWeight);
}

А вот это совсем не понятно
for (i = 0; i < m_FilterAperture; i++)
{
((COMPLEX_FLOAT32 *)m_Signal) = ((COMPLEX_FLOAT32 *)m_Signal)[i*m_Aperture/m_FilterAperture];
}
есть мнение что имелось ввиду
*((COMPLEX_FLOAT32 *)m_Signal) = ((COMPLEX_FLOAT32 *)m_Signal)*m_Aperture/m_FilterAperture
А вернее
for (i = 0; i < m_FilterAperture; i++)
{
((COMPLEX_FLOAT32 *)m_Signal).real = ((COMPLEX_FLOAT32 *)m_Signal).real*m_Aperture/m_FilterAperture;
((COMPLEX_FLOAT32 *)m_Signal).imag = ((COMPLEX_FLOAT32 *)m_Signal).imag*m_Aperture/m_FilterAperture;
}

505
17 февраля 2007 года
vAC
343 / / 28.02.2006
Не вижу никакой разницы между:
memset(&((COMPLEX_FLOAT32*)m_Signal)[Height], 0, (m_Aperture-Height)*sizeof(COMPLEX_FLOAT32));
и
memset(((COMPLEX_FLOAT32*)m_Signal)+Height, 0, (m_Aperture-Height)*sizeof(COMPLEX_FLOAT32));

Опять извиняюсь:
for (i = 0; i < m_FilterAperture; i++)
{
((COMPLEX_FLOAT32 *)m_Signal) = ((COMPLEX_FLOAT32 *)m_Signal)[i*m_Aperture/m_FilterAperture];
}
- это я пробовал разные варианты.

*((COMPLEX_FLOAT32 *)m_Signal) = ((COMPLEX_FLOAT32 *)m_Signal)*m_Aperture/m_FilterAperture
вообще работать не будет т.к. COMPLEX_FLOAT32 это не указатель.

И вообще-то не хотелось устраивать здесь уроки мастерства программирования, а хотелось бы по-сути...так какие мысли насчет преобразования такой простой последовательности?
563
17 февраля 2007 года
MrLinker
249 / / 17.09.2006
Цитата:
Точнее проблему сформулировать так:
как из частотной характеристики из N отсчетов получить характеристику состоящей M отсчетов, M<N, чтобы при обратном преобразовании небыло никаких затемнений...



Забавно..

Всё очень просто.
1. БПФ (N)
2. Играемся с коэф.
3. ОБПФ (M)

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

Т.е.
1. БПФ (N), тут получим N коэф в диапазоне 0...Fs/2
2. Играемся ...
3. ОБПФ (M), тут мы укорачиваем спектр и полагаем, что укороченный спектр опять занимает диапазон 0...Fs/2, что является бредовой идеей.

Очень туго я понимаю, что ты хочешь в итоге.

Может надо апроксимирировать спектр, сжимать в 16 раз и делать ОБПФ (M = N/16). Получим укороченный сигнал, по спектральному составу схожий с исходником.

Но в чем прикол??

А.. апроксимация не нужна, т.к. укорачивание происходит в Х раз, где Х - кратно 2. Т.е. делаем такое хитрое укорачивание исходного спектра, оставляя только каждую X*i гармонику.
Хотя...

505
17 февраля 2007 года
vAC
343 / / 28.02.2006
Я-то ничего не хочу получить, просто так было сделано до меня, а моя задача сделать быструю версию программы, и "темные" полоски мешают получить приемлемый результат. Говорю как все реализовано, сам к этому отношения не имел, поэтому ничего конкретного по поводу "зачем так, а не по другому" сказать не могу, сильно менять код не рекомендуется, а то еще результат будет неудовлетворительный и все шишки на меня посыпяться. Давайте тогда уж объясню в чем дело, не вдаваясь в подробности:

В обычной версии обрабатывается изображение целыми строками. У меня же строка разбивается на несколько частей и делается тоже самое над этими кусочками. И когда я их "склеиваю" получаются эти самые стыки в виде темных полос. В исходной версии затемненным оказывается только первый столбец, что не критично, а у меня он циклически повторяется в начале каждого маленького блока. Вот я и думаю, как мне избавиться от этого с минимальными изменениями в коде.
391
18 февраля 2007 года
Archie
562 / / 03.02.2005
Ну как вы не поймете, что фурье-преобразование корретно определено только для периодическиго сигнала. Если сигнал не периодический, то он дополняется до периодического путем копирования. Т.о. на местах стыков фрагментов возникают всплески (читай, ВЧ компоненты спектра значительной апмлитуды), которые потом зарезаются фильтром. В результате после восстановления сигнала на местах стыков - провалы (ВЧ частоты-то были вырезаны). Чтобы преодолеть проблему сию используют перекрытие фрагментов анализируемого сигнала (см. реализацию МП3 или любого другого кодека в частотном базисе).
505
18 февраля 2007 года
vAC
343 / / 28.02.2006
Цитата: Archie
Ну как вы не поймете, что фурье-преобразование корретно определено только для периодическиго сигнала. Если сигнал не периодический, то он дополняется до периодического путем копирования. Т.о. на местах стыков фрагментов возникают всплески (читай, ВЧ компоненты спектра значительной апмлитуды), которые потом зарезаются фильтром. В результате после восстановления сигнала на местах стыков - провалы (ВЧ частоты-то были вырезаны). Чтобы преодолеть проблему сию используют перекрытие фрагментов анализируемого сигнала (см. реализацию МП3 или любого другого кодека в частотном базисе).



Вот именно такое решение-то я и реализовал, перекрываю блоки. Другого ничего в голову-то и не пришло.

391
18 февраля 2007 года
Archie
562 / / 03.02.2005
Да другого ничего и нет, только окна и перекрытие фрагментов :)
26K
25 февраля 2007 года
Marty
5 / / 25.02.2007
Цитата: Archie
Да другого ничего и нет, только окна и перекрытие фрагментов :)



Согласен с Archie. Использование перекрытий - самый приемлемый вариант в сложившейся ситуации. Проще всего использовать половинное перекрытие с синусоидальным окном.

26K
25 февраля 2007 года
Marty
5 / / 25.02.2007
Цитата: MrLinker

Может надо апроксимирировать спектр, сжимать в 16 раз и делать ОБПФ (M = N/16). Получим укороченный сигнал, по спектральному составу схожий с исходником.

Но в чем прикол??



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

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

505
25 февраля 2007 года
vAC
343 / / 28.02.2006
Дело в том, что в исходном изображении присутствуют "волны" (скажем с периодом в 16 отсчетов), которые и надо ликвидировать. А затемнение вообщим-то из-за окна Кайзера и в тоже время нет. Дело в том, что если не применять окно, то появляются "засветления" :) Причем их протяженность намного больше, чем "затемнения" при взвешивании.
26K
26 февраля 2007 года
Marty
5 / / 25.02.2007
Цитата: vAC
Дело в том, что в исходном изображении присутствуют "волны" (скажем с периодом в 16 отсчетов), которые и надо ликвидировать. А затемнение вообщим-то из-за окна Кайзера и в тоже время нет. Дело в том, что если не применять окно, то появляются "засветления" :) Причем их протяженность намного больше, чем "затемнения" при взвешивании.



Насколько мне известно из курсов ЦОС использование окна Кайзера неудачно при реализации преобразования с перекрытиями из-за неоднозначного восстановления.
"Засветления" при отсутствии какого-либо окна у Вас появляются на интервале равным величине перекрытия (или удвоенной, не сообразить с ходу)? Тогда оно возникает скорее всего из-за неправильного взвешивания в области перекрытия - происходит сложение отсчетов и аплитуды увеличиваются. Когда при таком подходе сначала используется окно Кайзера, то "увеличение амплитуды" скрадывается и может даже получаться обратный эффект - затемнение.
Попробуйте все-таки использовать половинное перекрытие и синусоидальное окно, значение весов которого можно получить с пом. упрощенной функции:

void SinWindowFloat(float* win, int N)
{
int i;
double pi = 3.14159265358979;

for( i = 0; i < N; i++ )
{
win = (float)sin(pi/N*(i+0.5));
}
}

505
26 февраля 2007 года
vAC
343 / / 28.02.2006
Да с перекрытиями все нормально работает :)
26K
01 марта 2007 года
Marty
5 / / 25.02.2007
Цитата: vAC
Да с перекрытиями все нормально работает :)



Т.е. проблема решена и затемнений нет?

Реклама на сайте | Обмен ссылками | Ссылки | Экспорт (RSS) | Контакты
Добавить статью | Добавить исходник | Добавить хостинг-провайдера | Добавить сайт в каталог