Проблема с функцией создания октавного фильтра octaveFilter
Программа спектрального анализа данных лётных испытаний по шуму на местности. В программе используется функция octaveFilter, но почему-то она не работает в цикле так, как должна
Здравствуйте. Пишу сейчас диплом по программе спектрального анализа, одно из заданий - ускорение работы программы. 1/3-октавный анализ обычно длится около 10 минут, а 1/12-октавный 35+ минут, что при большом количестве записей очень затратно.
С помощью Profiler проанализировал код, 90-95% общего времени работы программы занимает функция фильтрации (приложу скриншот). Она выполняется в двойном цикле, где внешний цикл - это цикл по времени (вообще внешних цикла два - время разделено на две части, но в целом разница между двумя циклами незначительная), а внутренний - по частотам. Так как основное время занимает функция фильтрации, решил использовать параллельные вычисления, по-моему, другим способом ускорить вычисления и не получится.
Для параллельных вычислений использовал parfor. После исправления всех проблем программа запустилась и довольно быстро отработала (для 1/3-октавного анализа получилось 4 минуты вместо 10 в оригинале), но выявилась разница в массиве уровней звукового давления, если сравнить её с массивом оригинала. Где-то эта разница небольшая, во втором знаке, но где-то она достигает нескольких дБ, что неприемлемо (приложу скриншот с массивом-разницей между оригинальным массивом SPL и с parfor).
Первая строчка (первая итерация внешнего цикла по времени) считается правильно, а вот дальше начинаются проблемы. Видимо, происходит какое-то "подмешивание" данных других итераций, что странно, учитывая то, что при parfor они вроде как должны выполняться одновременно. После многих отладок выяснил, что проблема происходит как раз в функции фильтрации - на вход данные поступают одни и те же, но на выходе данные отличаются.
Обычно ещё в начале программы создаётся банк фильтров через цикл с использованием функции octaveFilter (*1), затем во внутреннем цикле выбирается фильтр конкретной частоты и происходит фильтрация (*2).
*1 Создание банка фильтров
for i=1:Nfc
octFiltBank{i} = octaveFilter('FilterOrder', N, 'CenterFrequency', F01(i), 'Bandwidth', BW, 'SampleRate', Fs, 'Oversample', true); %#ok
end*2 Внешний (по времени) и внутренний (по частотам) циклы без parfor (можно просто заменить for на parfor)
for ii=(2*W+1):l_spl-1 % Работа с оставшимися N-W секундами записи
xw = Dat_dec((Tab_sec_2(ii-2*W,1)+Tab_sec_2(ii-2*W,2)):(Tab_sec_2(ii,1)+Tab_sec_2(ii,2)+Max_lag-1),2);% Заготовка места под массив для фильтрации с запасом на максимальный фазовый сдвиг
for i=1:Nfc
octFilt = octFiltBank{i}; % Выборка фильтра из банка
yw = octFilt(xw); % Фильтрация
dMax_lag=Max_lag-LAG_121(i); % определение разницы между максимальным фазовым сдвигом и фактическим
yw(1:end-Max_lag)=yw(LAG_121(i)+1:end-dMax_lag); % устранение фазового сдвига
ywc=yw(1:end-Max_lag); % Выборка окна для экспоненциального взвешивания
ksi=((1:size(ywc,1)).').*(1/Fs); % Шкала времени для эксоненциальной взвешивающей функции
P2=ywc.^2;
EE=exp((ksi-size(ywc,1).*(1/Fs))./TAU);
SPL(ii,i)=20.*log10((((1./TAU).*(1/Fs).*sum((P2).*(EE))).^0.5)./0.000020); % Расчёт SPL по МЭК 61672-1
end
endКак-то раз решил попробовать создавать фильтр прямо во внутреннем цикле перед фильтрацией *3 (да, медленнее, но вдруг исправило бы проблему).
*3
for i=1:Nfc
octFilt = octaveFilter('FilterOrder', N, 'CenterFrequency', F01(i), 'Bandwidth', BW, 'SampleRate', Fs, 'Oversample', true); % создание фильтра
yw = octFilt(xw); % Фильтрация
dMax_lag=Max_lag-LAG_121(i); % определение разницы между максимальным фазовым сдвигом и фактическим
yw(1:end-Max_lag)=yw(LAG_121(i)+1:end-dMax_lag); % устранение фазового сдвига
ywc=yw(1:end-Max_lag); % Выборка окна для экспоненциального взвешивания
ksi=((1:size(ywc,1)).').*(1/Fs); % Шкала времени для эксоненциальной взвешивающей функции
P2=ywc.^2;
EE=exp((ksi-size(ywc,1).*(1/Fs))./TAU);
SPL(ii,i)=20.*log10((((1./TAU).*(1/Fs).*sum((P2).*(EE))).^0.5)./0.000020); % Расчёт SPL по МЭК 61672-1
endОказалось, что даже если использовать такой метод в оригинальной программе (без parfor), то всё равно будет возникать разница в массиве SPL. Причём самое интересно то, что эта разница будет абсолютно такой же, как при использовании parfor. То есть проблема не в parfor, а в том, что функция octaveFilter почему-то отказывается нормально работать в подобном цикле. Как я уже говорил, проблема начинается именно в строчке с фильтрацией.
Сталкивался ли кто-то с подобной проблемой? Есть ли у неё решение, или может какая-то схожая функция создания октавных фильтров? Заранее спасибо
Комментарии
Если верить Вам, то по всей видимости проблема в строчках:
octFilt = octFiltBank{i}; % Выборка фильтра из банкаyw = octFilt(xw); % ФильтрацияМожет быть проблема в свойствах фильтра, которая возникает при передаче разных объектов в одну переменную. Попробуйте заменить эти строки на
yw = octFiltBank{i}(xw);чтобы использовать объекты напрямую.
1. Замена двух строчек на одну "yw = octFiltBank{i}(xw);" ничего не поменяла, проблема остаётся.
2. Эталонные значения получены исходной программой, которая уже давно используется. Вообще эта программа - программный аналог аппаратной реализации. Изначально эта работа выполняются прямо на измерительной аппаратуре, но по определённым причинам хотим перейти на программную реализацию. Для этого нужно увеличить скорость работы. Результат работы программы с очень хорошей точность совпадает с результатом обработки на измерительной системе
Еще бы хотелось понять как вы получили эталонные значения, с которыми сравниваете результаты.
Значит надо выполнять отладку. Другого не дано. На каждом шаге цикла делать прерывание и проверять значения всех массивов. Только так ошибку отловить можно.
Отладку в paffor так просто не выполнить, но подручными средствами получилось. Расхождения начинаются именно на моменте фильтрации. Входные данные одинаковы, выходные разные. Все остальные данные (которые не зависят от данных фильтрации) тоже одинаковы. Вывод - проблема в фильтре. Косвенно данный вывод подтверждает то, что абсолютно такое же расхождение получается при создании фильтра прямо во внутреннем массиве перед фильтрацией
парфор отлаживать можно обычным фором, а когда уж отладите добавлять пар.
ПС: я не располагаю в данный момент временем чтобы вникат, но вы фильтр делаете циклом по времени? Парфором? я верно понял?
Так в том и дело, что for работает нормально, а вот с parfor расхождение. В любом случае, отладка сделана, место возникновения расхождения найдено - момент фильтрации. Повторю то, что уже написал в теме: обычно банк фильтров создаётся в своём отдельном цикле. Пример есть в теме, "*1 Создание банка фильтров".
Далее, для расчёта уровней звукового давления есть двойной цикл. Внешний цикл по времени, внутренний по частотам, во внутреннем цикле происходит фильтрация. Пример *1 ниже.
*1
И вот если во внутреннем цикле просто поменять for на parfor (или можно это сделать во внешнем, расхождение изменится - даже увеличится), то происходит расхождение расчётов.
Косвенно подтверждает проблему в фильтре то, что если в обычный внутренний цикл с for (который работает отлично) вместо выбора фильтра из банка (как в *1) добавить создание фильтра прямо перед фильтрацией (ниже такой пример *2 с внутренним циклом), то расхождение будет абсолютно таким же, как если цикл будет обычным (с выбором фильтра), но с parfor (пример *1, но внутренний цикл с parfor вместо for).
*2
octFilt , создаваеый функцией octaveFilter - это объект?
где он уничтожается?
вы уверены что в каждой новой итерации в объекте не остаются к.л. данные от "старого" фильтра?
Конкретно в данном примере да, это объект, создаваемый функцией. Зачистка octFilt ничего не даёт.
Проблема остаётся, расхождение такое же. Конкретно создание фильтра прямо в цикле не нужно, оно к тому же ещё и дополнительное время занимает. Просто я думаю, что природа возникновения ошибки одна и у создания фильтра в цикле, и у простого использования parfor, так как они дают абсолютно одинаковое расхождение относительно оригинала
- отладить вариант с одним фильтром, а затем уж добавлять банк?
- парфор тяжко настроить. вы уверены что переменные, создаваемые в одной итерации не связаны с другими?
Вот вариант, где банк фильтров создаётся как обычно, но затем из него всегда берётся только фильтр №7.
Приложу файл-скриншот разницы между SPL с pafror и обычным for при использовании только фильтра №7 (как в коде выше). Как видно, разница в простом изменении for на parfor есть даже при использовании всегда одного фильтра. Если задать этот один фильтр вне цикла, то разница всё равно будет, но другая
Тут видимо без обращения в техподдержку Mathwoks не обойтись.
объект octFilt может быть временной переменной?
На сколько я понимаю, octFilt и является временной, так как используется исключительно в цикле. За ссылку спасибо, попробую с этой информацией поработать.
Я недавно пытался решить проблему parfor и попробовал вместо создания банка фильтров в цикле сделать много отдельных переменных с разными фильтрами. Это происходит до циклов с рассчётом SPL.
И затем во внутреннем цикле через if (можно и case сделать) выбирать рабочим нужный фильтр.
Да, громоздко, но какой-то неведомой для меня причине такое исполнение даёт практически двукратное ускорение расчёта - вместо ~2 минут получается около 55 секунд. Это практически столько же, сколько и с parfor (а если ещё и parfor с этим использовать, то вообще 20 секунд получается, но результат такой же ошибочный, как обычно). Причём данные SPL полностью сходятся с оригиналом
Почему так происходит? По идее же Матлаб отлично должен работать с матрицами/векторами, а тут и дополнительное время на if, и в целом строчка выбора фильтра в Profile давала минимальный вклад по времени
Тут дело еще в том, что парфор необходимо использовать:
- только при крайней необходимости,
- матлаб изначально заточен для работы с матрицами. В нем по сути нет переменных, это матрицы из одного элемента. И если использовать правильные приемы работы с матрицами (присущие матлабу) и никаких циклов, то матлаб сам распараллеливает как надо.
Я когда-то давно долго пытался ускорить решение дифура. Мне надо было решить до**га раз один и тот же дифур с разными параметрами. Несколько тысяч раз. Вот я и парился через парфор. И в принципе не без результатно. НО!!!
Потом уже через пару лет, когда уже не актуальна задача была, сел, на бумажке продумал и расписал все, и все параметры в дифур передал один раз в виде матрицы. И оно мне выдало то же самое решение но вместо 30-40 минут, секунд за 20. (компьютер был тот же, и матлаб тот же) Т.е. матлаб сам хорошо распараллеливает внутри некоторых функций, главное правильно юзать функции. Но наверно не все функции так умеют. Надо подробно мануал читать.
UPD: https://docs.exponenta.ru/R2019a/parallel-computing/objects-and-handles-in-parfor-loops.html