Code covered by the BSD License  

Highlights from
Introduction to Simulink for Signal Processing

image thumbnail

Introduction to Simulink for Signal Processing

by

 

17 Jan 2008 (Updated )

Demo files that goes with the seminar with the same title.

myLMS_fixpt.m
%% Creating and Implementing LMS Algorithm to Filter out Cockpit Noise

%Turn on Logging
fipref('LoggingMode', 'on', 'DataTypeOverride', 'ForceOff');
%fipref('LoggingMode', 'on', 'DataTypeOverride', 'truedoubles');

% Define datawidths and fractions as varialbes in the data dictionary
iowidth = 16;
iofraclength = 15; 
prodsumwidth = 32;
prodsumfraclength = 29; 

% FIMATH describes the arithmetic operators work with FI Objects
F = fimath(...
        'RoundMode', 'nearest', ...
        'OverflowMode', 'saturate', ...
        'ProductMode','specifyprecision', ...
        'SumMode','specifyprecision', ...
        'ProductWordLength', prodsumwidth, 'SumWordLength', prodsumwidth, ...
        'ProductFractionLength',prodsumfraclength,'SumFractionLength',...
        prodsumfraclength,'RoundMode', 'floor');
    
 Tio = numerictype 


%% Load Pre-saved Test Voice Data

% Load Test Audio Signal to be used as clean input and play this sample
% through the PC sound card
load mtlb
sound(mtlb, Fs);
pause(2);

mtlb_fix = fi(mtlb,1,16, F); %Use pre-saved audio

%Display Log data for this variable
resetlog(mtlb_fix)
%[maxlog(mtlb_fix), minlog(mtlb_fix), noverflows(mtlb_fix), nunderflows(mtlb_fix)]


sound(double(mtlb_fix), Fs);

%Error between the fixed point and floating point values 
norm(double(mtlb_fix)-mtlb)

%% Create Random noise and filter it with an 'unknown' filtering process

% Observation noise signal corresponding to Engine noise
noise = randn(size(mtlb));
noise_fix = fi(noise,1,16,F); % Use the pre-defined FIMATH properties
range(noise_fix)

% 'Unknown' FIR system to be identified
b  = fir1(39,[.3 .5]);

% Convert this filter to DFILT Object to make analysis and fixed-point conversion easier
Hb = dfilt.dffir(b)
%Set Arithmetic as Fixed-Point
Hb.Arithmetic = 'fixed';
get(Hb)
info(Hb)

noise_filtered  = filter(b,1,noise);
noise_filtered_fix = filter(Hb, noise_fix); % Convert to Fixed Point

% Report the Logging Data for the filter
R = qreport(Hb)

%fipref('LoggingMode', 'on', 'DataTypeOverride', 'truedoubles');
mtlb_noisy_fix  = fi(mtlb(1:length(noise_filtered)),F) + fi(noise_filtered_fix,F);
dispLog('mtlb_noisy_fix');
%fipref('LoggingMode', 'on', 'DataTypeOverride', 'ForceOff');

% norm(-mtlb_noisy + double(mtlb_noisy_fix)) 

%% Plot the noisy spectrogram
mtlb_noisy_float = double(mtlb_noisy_fix);
subplot(1,2,1), spectrogram(mtlb_noisy_float,  128,120,256,'yaxis'); title('Noisy Audio')

%Play noisy signal
sound(double(mtlb_noisy_fix), Fs)
pause(2)

%% Create and Implement LMS Adaptive Filter
% Create simple LMS/nLMS/RLS Filter to identify filtering process and
% remove the filtered noise from desired signal

% Define Adaptive Filter Parameters
filterLength = 40;
weights = fi(zeros(1,filterLength),F);
step_size = fi(0.009,F);


% Initialize Filter's Operational inputs
output = fi(zeros(1, length(mtlb_noisy_fix)),1,16,13,F);
err = fi(zeros(1,length(mtlb_noisy_fix)),1,16,13,F);
input = zeros(1,filterLength);%,1,16,F);

dbstop if warning fi:overflow

% For Loop to run through the data and filter out noise
for n = 1: length(mtlb_noisy_fix)

    %Get input vector to filter
    for k= 1:filterLength
        if ((n-k)>0)
            input(k) = noise_fix(n-k+1);
        end
    end

    output(n) = weights * input';  %            <Output of Adaptive Filter
    %if n == 5, dispLog('output'), end
    %dispLog('output')
    
    err(n)  = mtlb_noisy_fix(n) - output(n); %          <Error Computation
    %if n == 5, dispLog('err'), end
    %dispLog('err')
 
    weights = weights + step_size * err(n) * input; %           <Weights  Updation 
    %if n == 5, dispLog('weights'), end
    %dispLog('weights')
    
    %plot(weights); % Optional command to plot weights to visualize how
    %they change
end

%fipref('LoggingMode', 'on', 'DataTypeOverride', 'ForceOff');

dispLog('output')
dispLog('err')
dispLog('weights')

%% Plot the Clean Spectrogram
err_float = double(err);
subplot(1,2,2), spectrogram(err_float, 128,120,256,'yaxis');title('Filtered Audio')

%% Play the filtered audio 

% The cleaned up audio signal is in this case the error signal
sound(err_float)

%%

Contact us