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

1     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2     % MATLAB Design: Adaptive Noise Canceler algorithm using Least Mean Square 
3     % (LMS) filter implemented in MATLAB
4     %
5     % Key Design pattern covered in this example: 
6     % (1) Use of function calls
7     % (2) Function inlining vs instantiation knobs available in the coder
8     % (3) Use of system objects in the testbench to stream test vectors into the design
9     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
10    
11    %#codegen
12    
13    function [filtered_signal, y, fc] = mlhdlc_lms_fcn(input, desired, step_size, reset_weights)
14    %
15    % 'input'  : The signal from Exterior Mic which records the ambient noise.
16    % 'desired': The signal from Pilot's Mic which includes 
17    %            original music signal and the noise signal
18    % 'err_sig': The difference between the 'desired' and the filtered 'input'
19    %            It represents the estimated music signal (output of this block)
20    % 
21    % The LMS filter is trying to retrieve the original music signal ('err_sig') 
22    % from Pilot's Mic by filtering the Exterior Mic's signal and using it to 
23    % cancel the noise in Pilot's Mic. The coefficients/weights of the filter 
24    % are updated(adapted) in real-time based on 'input' and 'err_sig'.
25    
26    % register filter coefficients
27    persistent filter_coeff;
28    if isempty(filter_coeff)
29        filter_coeff = zeros(1, 40);
30    end
31    
32    % Variable Filter
33    % call 'tapped_delay_fcn' function on path to create 40-step tapped delay
34    delayed_signal = mtapped_delay_fcn(input);
35    
36    % apply filter coefficients 
37    weight_applied = delayed_signal .* filter_coeff;
38    
39    % call 'treesum' function on matlab path to sum up the results
40    filtered_signal = mtreesum_fcn(weight_applied);
41    
42    % Output estimated Original Signal
43    td = desired;
44    tf = filtered_signal;
45    
46    esig = td - tf;
47    y = esig;
48    
49    % Update Weights
50    % call 'update_weight_fcn' function on matlab path to 
51    % calculate the new weights
52    updated_weight = update_weight_fcn(step_size, esig, delayed_signal, ...
53                                       filter_coeff, reset_weights);
54    
55    % update filter coefficients register
56    filter_coeff = updated_weight;
57    
58    fc = filter_coeff;
59    
60    
61    function y = mtreesum_fcn(u)
62    %Implement the 'sum' function without a for-loop
63    %  y = sum(u);
64    
65    %  The loop based implementation of 'sum' function is not ideal for 
66    %  HDL generation and results in a longer critical path. 
67    %  A tree is more efficient as it results in
68    %  delay of log2(N) instead of a delay of N delay
69    
70    %  This implementation shows how to explicitly implement the vector sum in 
71    %  a tree shape to enable hardware optimizations.
72    
73    %  The ideal way to code this generically for any length of 'u' is to use 
74    %  recursion but it is not currently supported by MATLAB Coder
75    
76    
77    % NOTE: To instruct MATLAB Coder to compile an external function, 
78    % add the following compilation directive or pragma to the function code
79    %#codegen
80    
81    % This implementation is hardwired for a 40tap filter.
82    
83    level1 = vsum(u);
84    level2 = vsum(level1);
85    level3 = vsum(level2);
86    level4 = vsum(level3);
87    level5 = vsum(level4);
88    level6 = vsum(level5);
89    y = level6;
90    
91    
92    
93    function output = vsum(input)
94    
95    coder.inline('always');
96    
97    vt = input(1:2:end);
98        
99    for i = int32(1:numel(input)/2)
100       k = int32(i*2);
101       vt(i) = vt(i) + input(k);
102   end
103   
104   output = vt;
105   
106   function tap_delay = mtapped_delay_fcn(input)
107   % The Tapped Delay function delays its input by the specified number 
108   % of sample periods, and outputs all the delayed versions in a vector
109   % form. The output includes current input
110   
111   % NOTE: To instruct MATLAB Coder to compile an external function, 
112   % add the following compilation directive or pragma to the function code
113   %#codegen
114   
115   persistent u_d;
116   if isempty(u_d)
117       u_d = zeros(1,40);
118   end
119   
120   
121   u_d = [u_d(2:40), input];
122   
123   tap_delay = u_d;
124   
125   
126   function weights = update_weight_fcn(step_size, err_sig, delayed_signal, filter_coeff, reset_weights)
127   % This function updates the adaptive filter weights based on LMS algorithm
128   
129   %   Copyright 2007-2015 The MathWorks, Inc.
130   
131   % NOTE: To instruct MATLAB Coder to compile an external function, 
132   % add the following compilation directive or pragma to the function code
133   %#codegen
134   
135   step_sig = step_size .* err_sig;
136   correction_factor = delayed_signal .* step_sig;
137   updated_weight = correction_factor + filter_coeff;
138   
139   if reset_weights
140       weights = zeros(1,40);
141   else    
142       weights = updated_weight;
143   end

1     
2     %
3     
4     %   Copyright 2011-2015 The MathWorks, Inc.
5     
6     clear ('mlhdlc_lms_fcn');
7     % returns an adaptive FIR filter System object,
8     % HLMS, that computes the filtered output, filter error and the filter
9     % weights for a given input and desired signal using the Least Mean
10    % Squares (LMS) algorithm.
11    
12    stepSize = 0.01;
13    reset_weights =false;
14    
15    hfilt = dsp.FIRFilter;          % System to be identified
16    hfilt.Numerator = fir1(10, .25);
17    
18    rng('default'); % always default to known state  
19    x = randn(1000,1);                         % input signal
20    d = step(hfilt, x) + 0.01*randn(1000,1);    % desired signal
21    
22    hSrc = dsp.SignalSource(x);
23    hDesiredSrc = dsp.SignalSource(d);
24    
25    hOut = dsp.SignalSink;
26    hErr = dsp.SignalSink;
27    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
28    %Call to the design
29    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
30    while (~isDone(hSrc))
31        [y, e, w] = mlhdlc_lms_fcn(step(hSrc), step(hDesiredSrc), stepSize, reset_weights);
32        step(hOut, y);
33        step(hErr, e);
34    end
35    
36    figure('Name', [mfilename, '_plot']);
37    subplot(2,1,1), plot(1:1000, [d,hOut.Buffer,hErr.Buffer]);
38    title('System Identification of an FIR filter');
39    legend('Desired', 'Output', 'Error');
40    xlabel('time index'); ylabel('signal value');
41    subplot(2,1,2); stem([hfilt.Numerator.', w(end-10:end).']);
42    legend('Actual','Estimated');
43    xlabel('coefficient #'); ylabel('coefficient value');