Проблема с функцией создания октавных фильтров octaveFilter

@The_Doctor49 · Регистрация: 06.01.2023

Author24 — интернет-сервис помощи студентам

Здравствуйте. Пишу сейчас диплом по программе спектрального анализа, одно из заданий - ускорение работы программы. 1/3-октавный анализ обычно длится около 10 минут, а 1/12-октавный 35+ минут, что при большом количестве записей очень затратно.

С помощью Profiler проанализировал код, 90-95% общего времени работы программы занимает функция фильтрации (приложу скриншот). Она выполняется в двойном цикле, где внешний цикл - это цикл по времени (вообще внешних цикла два - время разделено на две части, но в целом разница между двумя циклами незначительная), а внутренний - по частотам. Так как основное время занимает функция фильтрации, решил использовать параллельные вычисления, по-моему, другим способом ускорить вычисления и не получится.

Для параллельных вычислений использовал parfor. После исправления всех проблем программа запустилась и довольно быстро отработала (для 1/3-октавного анализа получилось 4 минуты вместо 10 в оригинале), но выявилась разница в массиве уровней звукового давления, если сравнить её с массивом оригинала. Где-то эта разница небольшая, во втором знаке, но где-то она достигает нескольких дБ, что неприемлемо (приложу скриншот с массивом-разницей между оригинальным массивом SPL и с parfor).

Первая строчка (первая итерация внешнего цикла по времени) считается правильно, а вот дальше начинаются проблемы. Видимо, происходит какое-то "подмешивание" данных других итераций, что странно, учитывая то, что при parfor они вроде как должны выполняться одновременно. После многих отладок выяснил, что проблема происходит как раз в функции фильтрации - на вход данные поступают одни и те же, но на выходе данные отличаются.

Обычно ещё в начале программы создаётся банк фильтров через цикл с использованием функции octaveFilter (*1), затем во внутреннем цикле выбирается фильтр конкретной частоты и происходит фильтрация (*2).

Matlab M

*1 Создание банка фильтров
 
for i=1:Nfc
    octFiltBank{i} = octaveFilter('FilterOrder', N, 'CenterFrequency', F01(i), 'Bandwidth', BW, 'SampleRate', Fs, 'Oversample', true); 
end

Matlab M

*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 (да, медленнее, но вдруг исправило бы проблему).

Matlab M

*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 почему-то отказывается нормально работать в подобном цикле. Как я уже говорил, проблема начинается именно в строчке с фильтрацией.

Сталкивался ли кто-то с подобной проблемой? Есть ли у неё решение, или может какая-то схожая функция создания октавных фильтров? Заранее спасибо

@АлександрКом · 09.01.2023, 13:01

The_Doctor49, заметил, что если создать фильтр и два раза прогнать через него один и тот же массив данных, то результаты будут разными. Может, ошибка как-то с этим связана?

Размеры разных xw будут одинаковыми или нет? Если да, то на octFilt можно подать сразу подготовленную матрицу, где каждый столбец - это ваши xw.

@The_Doctor49 · 10.01.2023, 01:01 **[ТС]**

Размер может незначительно отличаться. Проблема не в xw, так как я проверял все данные, идущие на вход фильтра. Сравнивал даже xw до фильтра и после, разницы никакой

@ViktorArs · 10.01.2023, 08:16

А я например не понимаю почему размер может отличаться?

@АлександрКом · 10.01.2023, 11:03

The_Doctor49, ниже кусок кода:

Matlab M

x1 = 1 : 10;
x2 = 11 : 20;
 
octFilt = octaveFilter('FilterOrder', 8, 'CenterFrequency', 1e3, 'Bandwidth', '1/3 octave', 'SampleRate', 1e4, 'Oversample', true);
y1 = octFilt(x1);
y2 = octFilt(x2);
 
octFilt1 = octaveFilter('FilterOrder', 8, 'CenterFrequency', 1e3, 'Bandwidth', '1/3 octave', 'SampleRate', 1e4, 'Oversample', true);
y3 = octFilt1(x2);

octFilt и octFilt1 - одинаковые фильтры, y2 и y3 - результат фильтрации x2. Но, если вы посмотрите на y2 и y3, то увидите, что они получились разными. Я думаю, что связано это с тем, что прежде чем подать на octFilt вектор x2, на него был подан вектор x1, что, возможно, как-то изменило внутренние настройки фильтра.

@ViktorArs · 10.01.2023, 13:49

Это нюансы работы объекта.
В реальной жизни крайне сложно понять как возможно изменение параметров фильтра сигналом. Он же неадаптивный? Или я ошибаюсь?

@АлександрКом 1296 / 926 / 448 Регистрация: 21.10.2012 Сообщений: 2,601
	09.01.2023, 13:01	2
	The_Doctor49, заметил, что если создать фильтр и два раза прогнать через него один и тот же массив данных, то результаты будут разными. Может, ошибка как-то с этим связана? Размеры разных xw будут одинаковыми или нет? Если да, то на octFilt можно подать сразу подготовленную матрицу, где каждый столбец - это ваши xw. 0

@The_Doctor49 0 / 0 / 0 Регистрация: 06.01.2023 Сообщений: 2
	10.01.2023, 01:01 [ТС]	3
	Размер может незначительно отличаться. Проблема не в xw, так как я проверял все данные, идущие на вход фильтра. Сравнивал даже xw до фильтра и после, разницы никакой 0

@ViktorArs 104 / 85 / 22 Регистрация: 09.01.2018 Сообщений: 463
	10.01.2023, 08:16	4
	А я например не понимаю почему размер может отличаться? 0

@ViktorArs 104 / 85 / 22 Регистрация: 09.01.2018 Сообщений: 463
	10.01.2023, 13:49	6
	Это нюансы работы объекта. В реальной жизни крайне сложно понять как возможно изменение параметров фильтра сигналом. Он же неадаптивный? Или я ошибаюсь? 0

Опции темы