Фильтрация методом преобразования Фурье
В программе есть фильтр, схема работы следующая образом:
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, чтобы при обратном преобразовании небыло никаких затемнений...
memcpy(&pSignal[FilterAperture/2], &pSignal[Aperture-FilterAperture/2], FilterAperture/2*sizeof(COMPLEX_FLOAT32));
&-пропущен
4. Взвешивание "укороченного" спектра окном Кайзера
3. ОБПФ "укороченного" спектра
Порядок цифр и пунктов не совподает.
Чето не понятно как-то не к месту применяется окно Кайзера.
Насчет затемнения, ты когда выкидываешь часть спектора, с ним уходит энергия. Но не в этом суть, затемнение тогда было бы на всех отчетах.
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 раз меньшего, чем исходный). А насчет окна Кайзера это я сказать не могу зачем...так было сделано до меня, а моя задача состоит совсем в другом, сделать которую как раз и мешают эти полосы.
Вообщим при обработке тестовой последовательности следующим образом:
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
Здесь уже на первом "светлый", на втором "темный", а дальше более или менее.
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;
}
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 это не указатель.
И вообще-то не хотелось устраивать здесь уроки мастерства программирования, а хотелось бы по-сути...так какие мысли насчет преобразования такой простой последовательности?
как из частотной характеристики из 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 гармонику.
Хотя...
В обычной версии обрабатывается изображение целыми строками. У меня же строка разбивается на несколько частей и делается тоже самое над этими кусочками. И когда я их "склеиваю" получаются эти самые стыки в виде темных полос. В исходной версии затемненным оказывается только первый столбец, что не критично, а у меня он циклически повторяется в начале каждого маленького блока. Вот я и думаю, как мне избавиться от этого с минимальными изменениями в коде.
Вот именно такое решение-то я и реализовал, перекрываю блоки. Другого ничего в голову-то и не пришло.
Согласен с Archie. Использование перекрытий - самый приемлемый вариант в сложившейся ситуации. Проще всего использовать половинное перекрытие с синусоидальным окном.
Может надо апроксимирировать спектр, сжимать в 16 раз и делать ОБПФ (M = N/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));
}
}
Т.е. проблема решена и затемнений нет?