• Регистрация
The_Doctor49
The_Doctor49 0.00
н/д

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

 

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

 

Файлы

  • Разница SPL между оригинальной программой и с parfor.png
  • Анализ 1_3 окт.jpg

Теги

    06.01.2023

    Комментарии

    • alextip
      alextip+45.13
      6.01.2023 12:47

      Если верить Вам, то по всей видимости проблема в строчках:

      octFilt = octFiltBank{i}; % Выборка фильтра из банка
      yw = octFilt(xw); % Фильтрация 

      Может быть проблема в свойствах фильтра, которая возникает при передаче разных объектов в одну переменную. Попробуйте заменить эти строки на

      yw = octFiltBank{i}(xw);

      чтобы использовать объекты напрямую.

      • The_Doctor49
        The_Doctor490.00
        6.01.2023 16:07

        1. Замена двух строчек на одну "yw = octFiltBank{i}(xw);" ничего не поменяла, проблема остаётся.

        2. Эталонные значения получены исходной программой, которая уже давно используется. Вообще эта программа - программный аналог аппаратной реализации. Изначально эта работа выполняются прямо на измерительной  аппаратуре, но по определённым причинам хотим перейти на программную реализацию. Для этого нужно увеличить скорость работы. Результат работы программы с очень хорошей точность совпадает с результатом обработки на измерительной системе

      • alextip
        alextip+45.13
        6.01.2023 12:56

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

        • alextip
          alextip+45.13
          6.01.2023 19:07

          Значит надо выполнять отладку. Другого не дано. На каждом шаге цикла делать прерывание и проверять значения всех массивов. Только так ошибку отловить можно.

          • The_Doctor49
            The_Doctor490.00
            6.01.2023 19:27

            Отладку в paffor так просто не выполнить, но подручными средствами получилось. Расхождения начинаются именно на моменте фильтрации. Входные данные одинаковы, выходные разные. Все остальные данные (которые не зависят от данных фильтрации) тоже одинаковы. Вывод - проблема в фильтре. Косвенно данный вывод подтверждает то, что абсолютно такое же расхождение получается при создании фильтра прямо во внутреннем массиве перед фильтрацией

            • aBoomest
              aBoomest+929.89
              6.01.2023 20:20

              парфор отлаживать можно обычным фором, а когда уж отладите добавлять пар.

              ПС: я не располагаю в данный момент временем чтобы вникат, но вы фильтр делаете циклом по времени? Парфором? я верно понял?

              • The_Doctor49
                The_Doctor490.00
                6.01.2023 20:40

                Так в  том и дело, что for работает нормально, а вот с parfor расхождение. В любом случае, отладка сделана, место возникновения расхождения найдено - момент фильтрации. Повторю то, что уже написал в теме: обычно банк фильтров создаётся в своём отдельном цикле. Пример есть в теме, "*1 Создание банка фильтров".

                Далее, для расчёта уровней звукового давления есть двойной цикл. Внешний цикл по времени, внутренний по частотам, во внутреннем цикле происходит фильтрация. Пример *1 ниже.

                *1

                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

                И вот если во внутреннем цикле просто поменять for на parfor (или можно это сделать во внешнем, расхождение изменится - даже увеличится), то происходит расхождение расчётов. 

                Косвенно подтверждает проблему в фильтре то, что если в обычный внутренний цикл с for (который работает отлично) вместо выбора фильтра из банка (как в *1) добавить создание фильтра прямо перед фильтрацией (ниже такой пример *2 с внутренним циклом), то расхождение будет абсолютно таким же, как если цикл будет обычным (с выбором фильтра), но с parfor (пример *1, но внутренний цикл с parfor вместо for).


                *2

                 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
                • aBoomest
                  aBoomest+929.89
                  6.01.2023 21:23

                  octFilt , создаваеый функцией octaveFilter - это объект?
                  где он уничтожается?
                  вы уверены что в каждой новой итерации в объекте не остаются к.л. данные от "старого" фильтра?

                  • The_Doctor49
                    The_Doctor490.00
                    6.01.2023 21:37

                    Конкретно в данном примере да, это объект, создаваемый функцией. Зачистка octFilt ничего не даёт.

                            octFilt = octaveFilter('FilterOrder', N, 'CenterFrequency', F01(i), 'Bandwidth', BW, 'SampleRate', Fs, 'Oversample', true);
                            yw = octFilt(xw); 
                            clear octFilt

                    Проблема остаётся, расхождение такое же. Конкретно создание фильтра прямо в цикле не нужно, оно к тому же ещё и дополнительное время занимает. Просто я думаю, что природа возникновения ошибки одна и у создания фильтра в цикле, и у простого использования parfor, так как они дают абсолютно одинаковое расхождение относительно оригинала

          • aBoomest
            aBoomest+929.89
            6.01.2023 22:17

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

            • The_Doctor49
              The_Doctor490.00
              6.01.2023 23:16

              Вот вариант, где банк фильтров создаётся как обычно, но затем из него всегда берётся только фильтр №7.

              for i=1:Nfc
                      octFilt = octFiltBank{7};
                      %octFilt = octaveFilter('FilterOrder', N, 'CenterFrequency', F01(i), 'Bandwidth', BW, 'SampleRate', Fs, 'Oversample', true);
                      yw = octFilt(xw); 
                      dMax_lag=Max_lag-LAG_30(i);
                      yw(1:end-Max_lag)=yw(LAG_30(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);
                  end

              Приложу файл-скриншот разницы между SPL с pafror и обычным for при использовании только фильтра №7 (как в коде выше). Как видно, разница в простом изменении for на parfor есть даже при использовании всегда одного фильтра. Если задать этот один фильтр вне цикла, то разница всё равно будет, но другая

            • alextip
              alextip+45.13
              7.01.2023 17:26

              Тут видимо без обращения в техподдержку Mathwoks не обойтись.

              • aBoomest
                aBoomest+929.89
                8.01.2023 21:23

                объект octFilt может быть временной переменной?

                • The_Doctor49
                  The_Doctor490.00
                  10.01.2023 01:18

                  На сколько я понимаю, octFilt и является временной, так как используется исключительно в цикле. За ссылку спасибо, попробую с этой информацией поработать.



                  Я недавно пытался решить проблему parfor  и попробовал вместо создания банка фильтров в цикле сделать много отдельных переменных с разными фильтрами. Это происходит до циклов с рассчётом SPL.

                  octFilt1 = octaveFilter('FilterOrder', N, 'CenterFrequency', F01(1), 'Bandwidth', BW, 'SampleRate', Fs, 'Oversample', true); 
                  octFilt2 = octaveFilter('FilterOrder', N, 'CenterFrequency', F01(2), 'Bandwidth', BW, 'SampleRate', Fs, 'Oversample', true); 
                  octFilt3 = octaveFilter('FilterOrder', N, 'CenterFrequency', F01(3), 'Bandwidth', BW, 'SampleRate', Fs, 'Oversample', true); 
                  
                  ...
                  
                  octFilt30 = octaveFilter('FilterOrder', N, 'CenterFrequency', F01(30), 'Bandwidth', BW, 'SampleRate', Fs, 'Oversample', true); 

                  И затем во внутреннем цикле через if (можно и case сделать) выбирать рабочим нужный фильтр.

                  for i=1:Nfc
                          
                           if(i==1)
                              octFiltq = octFilt1;
                           end
                          
                           if(i==2)
                              octFiltq = octFilt2;
                           end
                  
                               ...
                  
                           if(i==30)
                               octFiltq = octFilt30;
                           end
                          
                           yw = octFiltq(xw);    
                           dMax_lag=Max_lag-LAG_30(i);
                           yw(1:end-Max_lag)=yw(LAG_30(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);
                  end

                  Да, громоздко, но какой-то неведомой для меня причине такое исполнение даёт практически двукратное ускорение расчёта - вместо ~2 минут получается около 55 секунд. Это практически столько же, сколько и с parfor (а если ещё и parfor с этим использовать, то вообще 20 секунд получается, но результат такой же ошибочный, как обычно). Причём данные SPL полностью сходятся с оригиналом

                  Почему так происходит? По идее же Матлаб отлично должен работать с матрицами/векторами, а тут и дополнительное время на if, и в целом строчка выбора фильтра в Profile давала минимальный вклад по времени

                  • aBoomest
                    aBoomest+929.89
                    10.01.2023 14:16

                    Тут дело еще в том, что парфор необходимо использовать:
                    - только при крайней необходимости,
                    - матлаб изначально заточен для работы с матрицами. В нем по сути нет переменных, это матрицы из одного элемента. И если использовать правильные приемы работы с матрицами (присущие матлабу) и никаких циклов, то матлаб сам распараллеливает как надо.
                    Я когда-то давно долго пытался ускорить решение дифура. Мне надо было решить до**га раз один и тот же дифур с разными параметрами. Несколько тысяч раз. Вот я и парился через парфор. И в принципе не без результатно. НО!!!
                    Потом уже через пару лет, когда уже не актуальна задача была, сел,  на бумажке продумал и расписал все, и все параметры в дифур передал один раз в виде матрицы. И оно мне выдало то же самое решение но вместо 30-40 минут, секунд за 20. (компьютер был тот же, и матлаб тот же) Т.е. матлаб сам хорошо распараллеливает внутри некоторых функций, главное правильно юзать функции. Но наверно не все функции так умеют. Надо подробно мануал читать.

                • aBoomest
                  aBoomest+929.89
                  9.01.2023 10:39

                  Ближайшие события