Exponenta Banner
exponenta event banner

Ускорение алгоритма обработки сигналов в MATLAB

Этот пример использует:

Примечание

Бенчмарки в этом примере были измерены на машине с четырьмя физическими ядрами.

Этот пример показывает, как ускорить алгоритм обработки сигналов в MATLAB® использование codegen (MATLAB Coder) и dspunfold функций. Можно сгенерировать исполняемый файл MATLAB (MEX-функция) из всей функции MATLAB или определенных частей функции MATLAB. Когда вы запускаете MEX-функцию вместо оригинального кода MATLAB, скорость симуляции может значительно увеличиться. Чтобы сгенерировать эквивалент MEX, алгоритм должен поддержать генерацию кода.

Использовать codegen (MATLAB Coder), необходимо установить MATLAB Coder™. Использовать dspunfoldнеобходимо установить MATLAB Coder и DSP System Toolbox™.

Использовать dspunfold в Windows и Linux необходимо использовать компилятор, поддерживающий интерфейс приложения Open Multi-Processing (OpenMP). См. «Поддерживаемые компиляторы».

Алгоритм конечной импульсной характеристики

Рассмотрим простой алгоритм конечной импульсной характеристики, чтобы ускориться. Скопируйте firfilter код функции в firfilter.m файл.

function [y,z1] = firfilter(b,x)
% Inputs:
%   b - 1xNTaps row vector of coefficients
%   x - A frame of  noisy input 

% States:
%   z, z1 - NTapsx1 column vector of states

% Output:
%   y - A frame of filtered output
 
persistent z;

if (isempty(z))
    z = zeros(length(b),1);
end
Lx = size(x,1);
y = zeros(size(x),'like',x);

z1 = z;
for m = 1:Lx
    % Load next input sample
    z1(1,:) = x(m,:);
    
    % Compute output
    y(m,:) = b*z1;
    
    % Update states
    z1(2:end,:) = z1(1:end-1,:);
    z = z1;
end

 firfilter функция принимает вектор коэффициентов фильтра , b, шумный входной сигнал, x, в качестве входов. Сгенерируйте коэффициенты фильтра, используя fir1 функция.

NTaps = 250;
Fp = 4e3/(44.1e3/2);
b = fir1(NTaps-1,Fp);

Фильтруйте поток сигнала шумной синусоиды при помощи firfilter функция. Синусоида имеет формат кадра 4000 выборок и частоту дискретизации 192 кГц. Сгенерируйте синусоиду, используя dsp.SineWave Системные object™. Шум является белым Гауссовым со средним значением 0 и отклонением 0,02. Назовите эту функцию firfilter_sim.  firfilter_sim функция вызывает firfilter функция на шумном входе.

function totVal = firfilter_sim(b)
% Create the signal source
Sig = dsp.SineWave('SamplesPerFrame',4000,'SampleRate',19200);
totVal = zeros(4000,500);
R  = 0.02;

clear firfilter;

% Iteration loop. Each iteration filters a frame of the noisy signal.
for i = 1 : 500
    trueVal = Sig();                        % Original sine wave 
    noisyVal = trueVal + sqrt(R)*randn;     % Noisy sine wave
    filteredVal = firfilter(b,noisyVal);    % Filtered sine wave
    totVal(:,i) = filteredVal;              % Store the entire sine wave
end

Управляемый firfilter_sim и измерьте скорость выполнения. Скорость выполнения изменяется в зависимости от машины.

tic;totVal = firfilter_sim(b);t1 = toc;
fprintf('Original Algorithm Simulation Time: %4.1f seconds\n',t1);
Original Algorithm Simulation Time:  7.8 seconds

Ускорение конечной импульсной характеристики с помощью codegen

Звонить codegen на firfilter, и сгенерировать его MEX эквивалент, firfilter_mex. Сгенерируйте и передайте коэффициенты фильтра и сигнал синусоиды как входы к firfilter функция.

Ntaps = 250;
Sig = dsp.SineWave('SamplesPerFrame',4000,'SampleRate',19200); % Create the Signal Source
R = 0.02;           
trueVal = Sig();                    % Original sine wave
noisyVal = trueVal + sqrt(R)*randn; % Noisy sine wave
Fp = 4e3/(44.1e3/2);
b = fir1(Ntaps-1,Fp);               % Filter coefficients

codegen firfilter -args {b,noisyVal}

В firfilter_sim функция, замена firfilter(b,noisyVal) вызов функции с firfilter_mex(b,noisyVal). Назовите эту функцию firfilter_codegen.

function totVal = firfilter_codegen(b)
% Create the signal source
Sig = dsp.SineWave('SamplesPerFrame',4000,'SampleRate',19200);
totVal = zeros(4000,500);
R  = 0.02;

clear firfilter_mex;

% Iteration loop. Each iteration filters a frame of the noisy signal.
for i = 1 : 500
    trueVal = Sig();                           % Original sine wave 
    noisyVal = trueVal + sqrt(R)*randn;        % Noisy sine wave
    filteredVal = firfilter_mex(b,noisyVal);   % Filtered sine wave
    totVal(:,i) = filteredVal;                 % Store the entire sine wave
end

Управляемый firfilter_codegen и измерьте скорость выполнения. Скорость выполнения изменяется в зависимости от машины.

tic;totValcodegen = firfilter_codegen(b);t2 = toc;
fprintf('Algorithm Simulation Time with codegen: %5f seconds\n',t2);
fprintf('Speedup factor with codegen: %5f\n',(t1/t2));
Algorithm Simulation Time with codegen: 0.923683 seconds
Speedup factor with codegen: 8.5531

Ускорение усиления приблизительно 8.5.

Ускорение конечной импульсной характеристики с помощью dspunfold

dspunfold функция генерирует многопоточный файл MEX, который может улучшить ускорение усиления еще больше.

dspunfold также генерирует однопоточный файл MEX и функцию самодиагностического анализатора. Многопоточный файл MEX использует многоядерную архитектуру процессора хоста-компьютера. Однопоточный файл MEX аналогичен файлу MEX, который codegen функция генерирует. Функция анализатора измеряет ускорение усиления многопоточного файла MEX по однопоточному файлу MEX.

Звонить dspunfold на firfilter и сгенерируйте его многопоточный эквивалент MEX, firfilter_mt. Определите длину состояния в выборках с помощью -f опция, которая может улучшить ускорение усиления дальше. -s auto инициирует автоматическое обнаружение длины состояния. Для получения дополнительной информации об использовании -f и -s опции, см. dspunfold.

dspunfold firfilter -args {b,noisyVal} -s auto -f [false,true]
State length: [autodetect] samples, Repetition: 1, Output latency: 8 frames, Threads: 4
Analyzing: firfilter.m
Creating single-threaded MEX file: firfilter_st.mexw64
Searching for minimal state length (this might take a while)
Checking stateless ... Insufficient
Checking 4000 samples ... Sufficient
Checking 2000 samples ... Sufficient
Checking 1000 samples ... Sufficient
Checking 500 samples ... Sufficient
Checking 250 samples ... Sufficient
Checking 125 samples ... Insufficient
Checking 187 samples ... Insufficient
Checking 218 samples ... Insufficient
Checking 234 samples ... Insufficient
Checking 242 samples ... Insufficient
Checking 246 samples ... Insufficient
Checking 248 samples ... Insufficient
Checking 249 samples ... Sufficient
Minimal state length is 249 samples
Creating multi-threaded MEX file: firfilter_mt.mexw64
Creating analyzer file: firfilter_analyzer.p

Инструмент автоматического обнаружения длины состояния обнаруживает точную длину состояния 259 выборки.

Вызовите функцию анализатора и измерьте ускорение усиления многопоточного файла MEX относительно однопоточного файла MEX. Предоставьте по меньшей мере две различные системы координат для каждого входного параметра анализатора. Системы координат складываются по первой размерности. Анализатор чередуется между этими системами координат, проверяя, что выходы совпадают. Отказ нескольких систем координат для каждого входа может снизить эффективность анализатора и привести к ложным положительным результатам верификации. 

firfilter_analyzer([b;0.5*b;0.6*b],[noisyVal;0.5*noisyVal;0.6*noisyVal]);

Analyzing multi-threaded MEX file firfilter_mt.mexw64. For best results, 
please refrain from interacting with the computer and stop other processes until the 
analyzer is done.
Latency = 8 frames
Speedup = 3.2x

firfilter_mt имеет коэффициент ускоренного усиления 3.2 при сравнении с однопоточным файлом MEX, firfilter_st. Чтобы увеличить скорость дальнейшего, увеличьте коэффициент повторения с помощью -r опция. Компромисс заключается в том, что задержка на выходе увеличивается. Используйте коэффициент повторения 3. Задайте точную длину состояния, чтобы уменьшить накладные расходы и увеличить скорость еще больше.

dspunfold firfilter -args {b,noisyVal} -s 249 -f [false,true] -r 3

State length: 249 samples, Repetition: 3, Output latency: 24 frames, Threads: 4
Analyzing: firfilter.m
Creating single-threaded MEX file: firfilter_st.mexw64
Creating multi-threaded MEX file: firfilter_mt.mexw64
Creating analyzer file: firfilter_analyzer.p

Вызовите функцию анализатора.

firfilter_analyzer([b;0.5*b;0.6*b],[noisyVal;0.5*noisyVal;0.6*noisyVal]);
Analyzing multi-threaded MEX file firfilter_mt.mexw64. For best results, 
please refrain from interacting with the computer and stop other processes 
until the analyzer is done.
Latency = 24 frames
Speedup = 3.8x

Коэффициент ускорения усиления 3.8или примерно в 32 раза больше скорости выполнения исходной симуляции.

Для этого конкретного алгоритма можно увидеть, что dspunfold генерирует высоко оптимизированный код, без необходимости записывать какой-либо код C or C++. Ускорение усиления масштабируется с количеством ядер на вашей хост-машине.

Функция filter конечной импульсной характеристики в этом примере является только иллюстративным алгоритмом, который легко понять. Можно применить этот рабочий процесс к любому из пользовательских алгоритмов. Если вы хотите использовать конечная импульсная характеристика, рекомендуется использовать dsp.FIRFilter Системный объект в DSP System Toolbox. Этот объект запускается намного быстрее, чем номера бенчмарков, представленные в этом примере, без необходимости генерации кода.

Алгоритм фильтра Калмана

Рассмотрим алгоритм фильтра Калмана, который оценивает сигнал синусоиды от шумного входа. Этот пример показывает эффективность фильтра Калмана с codegen и dspunfold.

Входной вход шумной синусоиды имеет формат кадра 4000 выборки и частоту дискретизации 192 кГц. Шум является белым Гауссовым со средним значением 0 и отклонением 0,02.

Функция filterNoisySignal вызывает kalmanfilter функция на шумном входе. 

type filterNoisySignal

function totVal = filterNoisySignal
% Create the signal source
Sig = dsp.SineWave('SamplesPerFrame',4000,'SampleRate',19200);
totVal = zeros(4000,500);
R  = 0.02;
clear kalmanfilter;
% Iteration loop to estimate the sine wave signal
for i = 1 : 500
    trueVal = Sig();                    % Actual values
    noisyVal = trueVal + sqrt(R)*randn; % Noisy measurements
    estVal = kalmanfilter(noisyVal);    % Sine wave estimated by Kalman filter
    totVal(:,i) = estVal;               % Store the entire sine wave
end
type kalmanfilter

function [estVal,estState] = kalmanfilter(noisyVal)
% This function tracks a noisy sinusoid signal using a Kalman filter
%
% State Transition Matrix
A = 1;
stateSpaceDim = size(A,1);

% Measurement Matrix
H = 1;
measurementSpaceDim = size(H,1);
numTsteps = size(noisyVal,1)/measurementSpaceDim;

% Containers to store prediction and estimates for all time steps
zEstContainer = noisyVal;
xEstContainer = zeros(size(noisyVal));

Q = 0.0001; % Process noise covariance
R = 0.02;   % Measurement noise covariance
persistent xhat P xPrior PPrior;

% Local copies of discrete states
if isempty(xhat)
    xhat = 5; % Initial state estimate
end

if isempty(P)
    P = 1; % Error covariance estimate
end

if isempty(xPrior)
    xPrior = 0;
end

if isempty(PPrior)
    PPrior = 0;
end

% Loop over all time steps
for n=1:numTsteps
    
    % Gather chunks for current time step
    zRowIndexChunk = (n-1)*measurementSpaceDim + (1:measurementSpaceDim);
    stateEstsRowIndexChunk = (n-1)*stateSpaceDim + (1:stateSpaceDim);
        
    % Prediction step
    xPrior = A * xhat;
    PPrior =  A * P * A' + Q;
    
    % Correction step. Compute Kalman gain.
    PpriorH   = PPrior * H';
    HPpriorHR  = H * PpriorH + R;
    KalmanGain = (HPpriorHR \ PpriorH')';
    KH  = KalmanGain * H;
    
    % States and error covariance are updated in the
    % correction step
    xhat = xPrior + KalmanGain * noisyVal(zRowIndexChunk,:) - ...
        KH * xPrior;
    P = PPrior - KH * PPrior;
    
    % Append estimates
    xEstContainer(stateEstsRowIndexChunk, :) = xhat;
    zEstContainer(zRowIndexChunk,:) = H*xhat;
   
end

% Populate the outputs
estVal = zEstContainer;
estState = xEstContainer;

end

Управляемый filterNoisySignal.m и измерьте скорость выполнения.

tic;totVal = filterNoisySignal;t1 = toc;
fprintf('Original Algorithm Simulation Time: %4.1f seconds\n',t1);
Original Algorithm Simulation Time: 21.7 seconds

Ускорение фильтра Калмана с помощью codegen

Вызовите codegen функция on kalmanfilter, и сгенерировать его MEX эквивалент, kalmanfilter_mex.

kalmanfilter функция требует шумной синусоиды в качестве входов.

Sig = dsp.SineWave('SamplesPerFrame',4000,'SampleRate',19200); % Create the Signal Source
R = 0.02;                           % Measurement noise covariance
trueVal = step(Sig);                % Actual values
noisyVal = trueVal + sqrt(R)*randn; % Noisy measurements
codegen -args {noisyVal} kalmanfilter.m

Заменить kalmanfilter(noisyVal) в filterNoisySignal функция с kalmanfilter_mex(noisyVal). Назовите эту функцию как filterNoisySignal_codegen

function totVal = filterNoisySignal_codegen
% Create the signal source
Sig = dsp.SineWave('SamplesPerFrame',4000,'SampleRate',19200);
totVal = zeros(4000,500);
R  = 0.02;
clear kalmanfilter_mex;
% Iteration loop to estimate the sine wave signal
for i = 1 : 500
    trueVal = Sig();                        % Actual values
    noisyVal = trueVal + sqrt(R)*randn;     % Noisy measurements
    estVal = kalmanfilter_mex(noisyVal);    % Sine wave estimated by Kalman filter
    totVal(:,i) = estVal; % Store the entire sine wave
end

Управляемый filterNoisySignal_codegen и измерьте скорость выполнения.

tic; totValcodegen = filterNoisySignal_codegen; t2 = toc;
fprintf('Algorithm Simulation Time with codegen: %5f seconds\n',t2);
fprintf('Speedup with codegen is %0.1f',t1/t2);
Algorithm Simulation Time with codegen: 0.095480 seconds
Speedup with codegen is 227.0

Алгоритм фильтра Калмана реализует несколько матричных умножений. codegen использует Базовые подпрограммы линейной алгебры (BLAS) библиотеки для выполнения этих умножений. Эти библиотеки генерируют высоко оптимизированный код, что дает ускоренный коэффициент усиления 227.

Ускорение фильтра Калмана с помощью dspunfold

Сгенерируйте многопоточный файл MEX с помощью dspunfold и сравните его эффективность с codegen.

Sig = dsp.SineWave('SamplesPerFrame',4000,'SampleRate',19200); 
% Create the signal source
R = 0.02;                           % Measurement noise covariance
trueVal = step(Sig);                % Actual values
noisyVal = trueVal + sqrt(R)*randn; % Noisy measurements
dspunfold kalmanfilter -args {noisyVal} -s auto
State length: [autodetect] frames, Repetition: 1, Output latency: 8 frames, Threads: 4
Analyzing: kalmanfilter.m
Creating single-threaded MEX file: kalmanfilter_st.mexw64
Searching for minimal state length (this might take a while)
Checking stateless ... Insufficient
Checking 1 frames ... Sufficient
Minimal state length is 1 frames
Creating multi-threaded MEX file: kalmanfilter_mt.mexw64
Creating analyzer file: kalmanfilter_analyzer.p

Вызовите функцию анализатора.

kalmanfilter_analyzer([noisyVal;0.01*noisyVal;0.05*noisyVal;0.1*noisyVal]);
Analyzing multi-threaded MEX file kalmanfilter_mt.mexw64. For best results, 
please refrain from interacting with the computer and stop other processes until 
the analyzer is done.
Latency = 8 frames
Speedup = 0.7x

kalmanfilter_mt имеет коэффициент скорости 0.7, что является потерей эффективности 30% при сравнении с однопоточным файлом MEX, kalmanfilter_st. Увеличьте коэффициент повторения до 3 чтобы увидеть, увеличивается ли эффективность. Кроме того, обнаружите длину состояния в выборках.

dspunfold kalmanfilter -args {noisyVal} -s auto -f true -r 3

State length: [autodetect] samples, Repetition: 3, Output latency: 24 frames, Threads: 4
Analyzing: kalmanfilter.m
Creating single-threaded MEX file: kalmanfilter_st.mexw64
Searching for minimal state length (this might take a while)
Checking stateless ... Insufficient
Checking 4000 samples ... Sufficient
Checking 2000 samples ... Sufficient
Checking 1000 samples ... Sufficient
Checking 500 samples ... Sufficient
Checking 250 samples ... Insufficient
Checking 375 samples ... Sufficient
Checking 312 samples ... Sufficient
Checking 281 samples ... Sufficient
Checking 265 samples ... Sufficient
Checking 257 samples ... Insufficient
Checking 261 samples ... Sufficient
Checking 259 samples ... Sufficient
Checking 258 samples ... Insufficient
Minimal state length is 259 samples
Creating multi-threaded MEX file: kalmanfilter_mt.mexw64
Creating analyzer file: kalmanfilter_analyzer.p

Вызовите функцию анализатора.

kalmanfilter_analyzer([noisyVal;0.01*noisyVal;0.05*noisyVal;0.1*noisyVal]);

Analyzing multi-threaded MEX file kalmanfilter_mt.mexw64. For best results, 
please refrain from interacting with the computer and stop other processes until the 
analyzer is done.
Latency = 24 frames
Speedup = 1.4x

dspunfold дает ускоренный коэффициент усиления 40% по сравнению с высоко оптимизированным однопоточным файлом MEX. Задайте точную длину состояния и увеличьте коэффициент повторения до 4. 

dspunfold kalmanfilter -args {noisyVal} -s 259 -f true -r 4

State length: 259 samples, Repetition: 4, Output latency: 32 frames, Threads: 4
Analyzing: kalmanfilter.m
Creating single-threaded MEX file: kalmanfilter_st.mexw64
Creating multi-threaded MEX file: kalmanfilter_mt.mexw64
Creating analyzer file: kalmanfilter_analyzer.p

Активируйте функцию анализатора, чтобы увидеть ускорение усиления.

kalmanfilter_analyzer([noisyVal;0.01*noisyVal;0.05*noisyVal;0.1*noisyVal]);

Analyzing multi-threaded MEX file kalmanfilter_mt.mexw64. For best results, please refrain 
from interacting with the computer and stop other processes until the analyzer is done.
Latency = 32 frames
Speedup = 1.5x

Коэффициент ускорения усиления 50% по сравнению с однопоточным файлом MEX.

Коэффициенты усиления эффективности codegen и dspunfold дайте зависимость от вашего алгоритма. codegen обеспечивает достаточное ускорение для некоторых конструкций MATLAB. dspunfold может обеспечить дополнительное увеличение эффективности с помощью ядер, доступных на вашей машине, чтобы распределить алгоритм через развертывание. Как показано в этом примере, количество ускорения, которое dspunfold обеспечивает зависит от конкретного алгоритма для ускорения. Использовать dspunfold в дополнение к codegen если ваш алгоритм хорошо подходит для распределения через развертывание, и если полученная в результате стоимость задержки соответствует ограничениям вашего приложения.

Некоторые конструкции MATLAB высоко оптимизированы с интерпретированным выполнением MATLAB. fft функция, например, запускается намного быстрее в интерпретированной симуляции, чем с генерацией кода.

Похожие темы

Документация DSP System Toolbox

Поддержка