I'm trying to convert Matlab code to C++ or C, but it has some errors. Can someone help me with this?

%function [e,w]=lmsBasic(M,x,d);
% Input arguments:
% M = learned filter length, dim 1x1
% x = input signal, dim Nx1
% d = desired signal, dim Nx1
% h = input coefficients  (coeff of filter we
%       are trying to learn. filter
%       that turns d into u)
%
% Output arguments:
% e = estimation error, dim Nx1
% w = final filter coefficients, dim Mx1
mus = [0.01 0.05 0.1 0.5 1 2];% 0.5 1];
for j = 1:length(mus)
    for i = 1:100
        %%INITIALIZE VALUES%%
        % generate input signal
        M=128; %buffer size (num filter weights)
        %=audioread('Airplane Sound Effect.m4a'); %input signal
        x=x/max(x); %sample rate
        fs=8000; %number of samples of the input signal
        N=length(x); %length of input signal
        muOG = mus(j);
        
        % generate known filter coefficients
        Pz=(0.5*[0:127]).^2;  %linear coefficients
        %Pz=randn(128,1); %random coefficients
        ylim = max(Pz)*1.20;
        ymin = min(Pz)-.2*max(Pz);
        % generate filtered input signal == desired signal
        d=conv(Pz,x); %input signal filtered by known filter Pz (primary path)
        
        
        %% LMS FOR MAIN ANC %%
        
        %initalize Wz filter values
        Wz=zeros(M,1);
        emean=zeros(N,1);
        %Make sure that x and d are column vectors
        x=x(:);
        d=d(:);
        %LMS
        
        for n=M:N
            xvec=x(n:-1:n-M+1); %input has to be in reverse orxer
            %adaptively update mu
            mu(n) = muOG;
            e(n)=d(n)-Wz'*xvec; %update error
            Wz=Wz+mu(n)*xvec*(e(n)); %update filter coefficient
            
            %visualize learned filter in realtime
            %         plot(Pz)
            %         hold on
            %         plot(Wz)
            %         axis([0 inf ymin ylim])
            %         title(sprintf('n=%f time=%fs error = %f mu=%f',n-M, (n-M)/fs, e(n), mu(n-M+1)))
            %         hold off
            %         legend('Input coefficients','Learned Coefficients')
            %         drawnow;
        end
        e=e(:);
        emean = (emean(:)+e);
    end
    emean=(emean)/i;
    if sum(isinf(emean))>0
        emean(~isinf(emean))=1e3;
        eall(j,:)=emean;
    elseif max(emean)>1e4
        for l = 1:length(emean)
            if abs(emean(l)) > 1e3
                emean(l)=1e3;
            end
        end
        eall(j,:)=emean;
    else
        [eall(j,:),q]=(envelope(abs(emean),150,'peaks'));
    end
end
figure
for i = 1:length(mus)
    plot(10*log10(abs(eall(i,:))))
    hold on
end
title('Convergence Time in Cycles')
ylabel('Error (dB)');
xlabel('Cycles');
hleg=legend('0.01','0.05','0.1','0.5','1.0', '2.0');
htitle=get(hleg,'Title');
set(htitle,'String','mu');
pause(17);
soundsc(x,fs);

% %% PLOT RESULTS %%
% figure
% subplot(2,1,1)
% plot(e)
% title('Convergence Time in Cycles')
% ylabel('Amplitude');
% xlabel('Cycles');
% legend('Error');
% subplot(2,1,2)
% stem(Pz)
% hold on
% stem(Wz, 'r*')
% title('Input Coefficients vs Learned Coefficients')
% ylabel('Amplitude');
% xlabel('Numbering of filter tap');
% legend('Input Coefficients', 'learned coefficients')

[Walter Roberson]

What error are you encountering?

Are you converting by hand, or are you using MATLAB Coder?

[Sean Roh]

I'm using MATLAB Coder. Can you try to convert this code in C++

[Walter Roberson]

You need a function statement for MATLAB Coder.

MATLAB Coder cannot process graphics or sound.

[Sean Roh]

I know that MATLAB Coder cannot process graphics or sound, so I'm going to find a way to do that. But can you actually try to convert the code because I have been working on this few days and I can't figure out how to convert it. You can comment out the the audioread part and graphing part.

[Sean Roh]

This is the original code.

%function [e,w]=lmsBasic(M,x,d,h);
%
% Input arguments:
% M = learned filter length, dim 1x1
% x = input signal, dim Nx1
% d = desired signal, dim Nx1
% h = input coefficients  (coeff of filter we
%       are trying to learn. filter
%       that turns d into u)
%
% Output arguments:
% e = estimation error, dim Nx1
% w = final filter coefficients, dim Mx1
mus = [0.01 0.05 0.1 0.5 1 2];% 0.5 1];
for j = 1:length(mus)
    for i = 1:100
        %%INITIALIZE VALUES%%
        % generate input signal
        M=128; %buffer size (num filter weights)
        x=randn(20000,1); %input signal
        x=x/max(x); %sample rate
        fs=8000; %number of samples of the input signal
        N=length(x); %length of input signal
        muOG = mus(j);
        
        % generate known filter coefficients
        Pz=(0.5*[0:127]).^2;  %linear coefficients
        %Pz=randn(128,1); %random coefficients
        ylim = max(Pz)*1.20;
        ymin = min(Pz)-.2*max(Pz);
        % generate filtered input signal == desired signal
        d=conv(Pz,x); %input signal filtered by known filter Pz (primary path)
        
        
        %% LMS FOR MAIN ANC %%
        
        %initalize Wz filter values
        Wz=zeros(M,1);
        emean=zeros(N,1);
        %Make sure that x and d are column vectors
        x=x(:);
        d=d(:);
        %LMS
        
        for n=M:N
            xvec=x(n:-1:n-M+1); %input has to be in reverse orxer
            %adaptively update mu
            mu(n) = muOG;
            e(n)=d(n)-Wz'*xvec; %update error
            Wz=Wz+mu(n)*xvec*(e(n)); %update filter coefficient
            
            %visualize learned filter in realtime
            %         plot(Pz)
            %         hold on
            %         plot(Wz)
            %         axis([0 inf ymin ylim])
            %         title(sprintf('n=%f time=%fs error = %f mu=%f',n-M, (n-M)/fs, e(n), mu(n-M+1)))
            %         hold off
            %         legend('Input coefficients','Learned Coefficients')
            %         drawnow;
        end
        e=e(:);
        emean = (emean(:)+e);
    end
    emean=(emean)/i;
    if sum(isinf(emean))>0
        emean(~isinf(emean))=1e3;
        eall(j,:)=emean;
    elseif max(emean)>1e4
        for l = 1:length(emean)
            if abs(emean(l)) > 1e3
                emean(l)=1e3;
            end
        end
        eall(j,:)=emean;
    else
        [eall(j,:),q]=(envelope(abs(emean),150,'peaks'));
    end
end
figure
for i = 1:length(mus)
    plot(10*log10(abs(eall(i,:))))
    hold on
end
title('Convergence Time in Cycles')
ylabel('Error (dB)');
xlabel('Cycles');
hleg=legend('0.01','0.05','0.1','0.5','1.0', '2.0');
htitle=get(hleg,'Title');
set(htitle,'String','mu');

% %% PLOT RESULTS %%
% figure
% subplot(2,1,1)
% plot(e)
% title('Convergence Time in Cycles')
% ylabel('Amplitude');
% xlabel('Cycles');
% legend('Error');
% subplot(2,1,2)
% stem(Pz)
% hold on
% stem(Wz, 'r*')
% title('Input Coefficients vs Learned Coefficients')
% ylabel('Amplitude');
% xlabel('Numbering of filter tap');
% legend('Input Coefficients', 'learned coefficients')

Answer 1

[Walter Roberson]

My adjusted version of the code passes compilation, and produces graphs. However I can demonstrate that the algorithm is wrong. And besides that, it produces NaN values.

function lmsBasic()
    %
    % Input arguments:
    % M = learned filter length, dim 1x1
    % x = input signal, dim Nx1
    % d = desired signal, dim Nx1
    % h = input coefficients  (coeff of filter we
    %       are trying to learn. filter
    %       that turns d into u)
    %
    % Output arguments:
    % e = estimation error, dim Nx1
    % w = final filter coefficients, dim Mx1
    
    M=128; %buffer size (num filter weights)
    
    mus = [0.01 0.05 0.1 0.5 1 2];% 0.5 1];
    fs=8000; %number of samples of the input signal
    Nrand = 20000;
    Niter = 100;
    
    Nmus = length(mus);
    
    mu = zeros(Nrand,1);
    e = zeros(Nrand,1);
    eall = zeros(Nmus, Nrand);
    
    for j = 1:Nmus
        emeans = zeros(Nrand, 1);
        for i = 1:Niter
            %%INITIALIZE VALUES%%
            % generate input signal
            x=randn(Nrand,1); %input signal
            x=x/max(x); %sample rate
            muOG = mus(j);
            
            % generate known filter coefficients
            Pz=(0.5*(0:127)).^2;  %linear coefficients
            %Pz=randn(128,1); %random coefficients
            ylim = max(Pz)*1.20;
            ymin = min(Pz)-.2*max(Pz);
            % generate filtered input signal == desired signal
            d=conv(Pz,x); %input signal filtered by known filter Pz (primary path)
            
            
            %% LMS FOR MAIN ANC %%
            
            %initalize Wz filter values
            Wz=zeros(M,1);
            emean=zeros(Nrand,1);
            %Make sure that x and d are column vectors
            x=x(:);
            d=d(:);
            %LMS
            
            for n=M:Nrand
                xvec=x(n:-1:n-M+1); %input has to be in reverse orxer
                %adaptively update mu
                mu(n) = muOG;
                e(n)=d(n)-Wz'*xvec; %update error
                Wz=Wz+mu(n)*xvec*(e(n)); %update filter coefficient
            end
            e=e(:);
            emeans = emean(:)+e;
        end
        emeans=(emeans)/Niter;
        if sum(isinf(emeans))>0
            emeans(~isinf(emeans))=1e3;
            eall(j,:)=emeans;
        elseif max(emeans)>1e4
            for l = 1:length(emeans)
                if abs(emeans(l)) > 1e3
                    emeans(l)=1e3;
                end
            end
            eall(j,:)=emeans;
        else
            [eall(j,:),q]= envelope(abs(emeans),150,'peaks');
        end
    end
    figure
    for i = 1:length(mus)
        plot(10*log10(abs(eall(i,:))))
        hold on
    end
    title('Convergence Time in Cycles')
    ylabel('Error (dB)');
    xlabel('Cycles');
    %{
hleg=legend('0.01','0.05','0.1','0.5','1.0', '2.0');
htitle=get(hleg,'Title');
set(htitle,'String','mu');
    %}

[Walter Roberson]

Referring to your original code:

    for i = 1:100

so i will be done 100 times.

        x=randn(20000,1); %input signal

x is 20000 long

        N=length(x); %length of input signal

so N will be 20000.

        for n=M:N

n can go to N and N is 20000

            e(n)=d(n)-Wz'*xvec; %update error

so e can be up to 20000 long.

        emean=zeros(N,1);

so emean is length 20000. emean is not updated inside for n so after for n it will still be zeros of length 20000. Note that emean is set to zero within for i so it will be reset to zero every iteration of for i

        e=e(:);

e is length 20000 as described above, and this statement just makes sure it becomes a column vector.

        emean = (emean(:)+e);

emean is, as described above, zeros(N,1) so when you add e to it, you get the same as e, 20000 x 1.

      end

That is the for i loop. And every iteration of i you reset emean to 0. The result is that when you come out of the for i loop, that emean only has data from the last iteration of the for i loop. There was no point in doing the rest of the iterations. Therefore the algorithm is wrong.

[Walter Roberson]

You will also find that when j = 4, mus 0.5, that about 9700 iterations into for n that the line

                e(n)=d(n)-Wz'*xvec; %update error

will lead to NaN. This arises because Wz'*xvec becomes inf as the Wz entries climb and climb until they are on the order of 1e307 and the sum of those times xvec (which ranges from about -0.6 to +0.6) overflows to infinity.

[Sean Roh]

Is this mex file?

[Walter Roberson]

I cannot test the mex file because I have not installed a compiler on my system. However, I can get it right through to the point where it tests the file at the MATLAB level, finds no run-time errors, and generates the mex file.

[Sean Roh]

Are you familiar with any real-time audio? Because I'm trying to make this real-time.

[Walter Roberson]

I do have some experience with real-time audio. However, as far as I know, I lost my extensive source-code a few decades ago; otherwise you might have been able to write a small parser to convert the RSX-11 assembly code into C, unless you happened to have a spare PDP11/23 on hand.

[Sean Roh]

Do you know anyone that can help me with real-time active noise cancellation?

[Walter Roberson]

https://www.mathworks.com/help/audio/ug/active-noise-control-using-a-filtered-x-lms-fir-adaptive-filter.html