Exponenta Banner

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

Этот пример использование:

Примечание

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

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

Чтобы использовать codegen, необходимо было установить 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. Синусоида имеет формат кадра 4 000 выборок и частоту дискретизации 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 или Код С++. Усиление ускорения масштабируется с количеством ядер на вашей хост-машине.

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

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

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

Шумный вход синусоиды имеет формат кадра 4 000 выборок и частоту дискретизации 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 на 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 пользуется библиотеками Basic Linear Algebra Subroutines (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
Поддержка