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

Алгоритм фильтра FIR

Рассмотрим простой алгоритм FIR-фильтра для ускорения. Копировать 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

Ускорение фильтра FIR с помощью 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.

Ускорение фильтра FIR с помощью 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 или C++. Увеличение скорости масштабируется с количеством ядер на хост-машине.

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

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

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