Programs/code

The following are privately generated scripts used for various test regimes, as well as plotting the results.

Matlab code 

fft_hz.m 

function [ frequency,amplitude, THD_percentage] = fft_hz( varargin ) 

%fft_hz( data,Fs ) can calculate multichannel sound samples 

%   frequecny in hz, amplitde in db FS, THD_percentage in % 

data=varargin{1}; 

Fs=varargin{2}; 

[num_of_samples,number_of_channel]=size(data); 

%preallocating 

y=zeros(num_of_samples,number_of_channel); 

amplitude=zeros(round(num_of_samples/2)+1,number_of_channel); 

total_harmonic_distortion=zeros(1,number_of_channel); 

frequency = Fs/2*linspace(0,1,round(num_of_samples/2)+1)'; 

for ii=1:number_of_channel 

% y(:,ii) = fft(data(:,ii),num_of_samples)/num_of_samples;% This can be used for better speed 

y(:,ii) = fft(data(:,ii))/num_of_samples; 

amplitude(:,ii) = 20*log10(abs(y(1:round(num_of_samples/2)+1,ii))); 

%amplitude(:,ii) = (abs(y(1:num_of_samples/2+1,ii))); 

total_harmonic_distortion(ii) = thd(data(:,ii),Fs,3); 

end 

THD_percentage=10.^(total_harmonic_distortion/20)*100; 

end 

audio_acqusition_from_daq.m 

%% Multichannel Audio Input and Output 

% this is the main code to for testing using Lindos microphone 

% the microphone input is differential input 

% for example the 'signal+' and 'signal-' output from one microphone  

% is going to the DAQ 'analog input channel_xx +' and 'analog  input channel_xx -' 

% for 2 microphones, it requires 2*2=4 analog input channels 

% the fifth analog input channel is to get amplifier output, which is used as a reference to 

% see if the output signal from the amplifier is clean enough. 

% the DAQ used analog output channel ao0 or ao1 to generate signal to the 

% amplifier 

% Clean up 

clc;clear all;close all; 

%% Specify sampling rate and duration 

Fs=200000;% specify sampling rate e.g.Fs=192000;%make sure that your soundcard supports it 

test_duration=0.2;% specify testing duration in [s] 

Fs_out=800e3; 

%% Choose input and output devices 

%%% Do this if it is the first time you run this code 

%%% 

%d = daq.getDevices % choose devices, important!!! make sure devices is selected correctly 

%% Create session for microphone input 

s = daq.createSession('ni') ; 

s.Rate=Fs; 

s.DurationInSeconds=test_duration; 

% add audio input channel 

%do this for reference microphone 

ch1=addAnalogInputChannel(s,'Dev1', [6,14,7,15,11], 'Voltage'); 

% note: 

% 6 and 14 is the + and - input channel 

% 7 and 15 is the + and - input channel 

% 11 is a single ended input, the signal is from the amplifier output 

number_of_input_channel=length(ch1); % generally 2 (left and right channel) 

for nn=1:number_of_input_channel 

    ch1(nn).Range=[-1,1]; % Range: -1.0 to +1.0 Volts; 

end 

%ch1(5).Range=[-5,5]; 

%% Create session for test signal output 

s2 = daq.createSession('ni'); 

s2.Rate=Fs_out; 

s.DurationInSeconds = test_duration; 

ch_out=addAnalogOutputChannel(s2,'Dev1', 0, 'Voltage'); % important, make sure this is selected correctly 

ch_out.Range=[-5,5]; 

number_of_output_channel=length(ch_out);% generally 2 (left and right channel) 

for nn=1:40 

    %% Generate test signal 

    outputDurationInSeconds = test_duration+0.1; 

    % Determine sample times. 

    test_time = (1/Fs_out)*[1:round(outputDurationInSeconds*Fs_out)]'; 

    % % Generate white noise test signal 

    % fc = (Fs_out/2)*rand(size(test_time)); A  = 0.2; 

    % test_signal = sin(2*pi*fc.*test_time); 

    % % Generate sinewave test tone. 

    fc = 500*nn; A = 3; 

    test_signal = A*sin(2*pi*fc*test_time); 

    % %    Generate linear sweep 

    %test_signal = 3*chirp(test_time,10,0.8*test_duration,20000);     % Start @ 15Hz, 

    %                               %   cross 20000Hz at t=0.8*test_duration [sec] 

    % figure() 

    % plot(test_time,test_signal) 

    %% outputting test_signal in background and getting input in foreground 

    queueOutputData(s2, repmat(test_signal, 1, number_of_output_channel));% repmat to match output data to channel 

    lh=addlistener(s2,'DataRequired', ... 

        @(src,event) src.queueOutputData(s2, repmat(y, 1, number_of_output_channel)));% add listener for output 

    startBackground(s2); 

    %pause(0.1); % matlab has an internal delay,pause to avoid getting useless signal input 

    [data,t] = startForeground(s); 

    delete (lh) % delete listener 

    %% bandpass filter 

    d = fdesign.bandpass('N,F3dB1,F3dB2',10,200,60e3,Fs); 

    Hd = design(d,'butter'); 

    % fvtool(Hd) 

    % data2=data; 

    % for ii=1:number_of_input_channel 

    % data(:,ii)=filter(Hd,data2(:,ii)); 

    % end 

    %% plot result 

    %do this if it is differential input. 

    %%%%%%%%%%%%%% 

    number_of_input_channel=3; 

    data_1=0.5*(data(:,1)-data(:,2)); 

    data_2=0.5*(data(:,3)-data(:,4)); 

    data=[data_1 data_2 data(:,5)]; 

    data2=data; 

    for ii=1:number_of_input_channel 

        data(:,ii)=filter(Hd,data2(:,ii)); 

    end 

    %%%%%%%%%%%%% 

    fprintf('Testing for frequency=%d [Hz]\n',fc) 

    title_str={'Reference 1','Reference 2','Amplifier output'}; 

    figure() 

    for ii=1:number_of_input_channel 

        subplot(2,2,ii); 

        plot(t(200:round(Fs/fc)*4+200), data(200:round(Fs/fc)*4+200,ii)); 

        %plot(t, data(:,ii)); 

        %         title([s.Channels(ii).Device.Model(1:3) '...' ' Ch ' s.Channels(ii).ID... 

        %             ]); 

        title(title_str(ii)); 

        %         title(['f=' num2str(fc)]) 

    end 

    % use fft_hz to get response and THD 

    [ frequency,amplitude, THD_percentage] = fft_hz( data,Fs); 

    %amplitude = sgolayfilt(amplitude,7,39);% curve fitting 

    % 20log10(A)+119 for lindos mm4 with 108db setting 

    % 20log10(A)+129 for lindos mm4 with 118db setting 

    % 20log10(A)+139 for lindos mm4 with 128db setting 

    sensitivity=[139,139,0]; 

    amplitude=repmat(sensitivity,length(amplitude),1)+amplitude; 

    [amp_desc, index] = sort(amplitude,'descend'); 

    figure(); 

    for ii=1:number_of_input_channel 

        %semilogx(frequency,amplitude(:,ii)) 

        if ii==1|| ii==2 

            frequency_plot_length=round(length(frequency)/4); 

            subplot(2,2,ii); 

        else 

            frequency_plot_length=length(frequency); 

            subplot(2,2,3:4); 

        end 

        plot(frequency(1:frequency_plot_length),amplitude(1:frequency_plot_length,ii),'k','LineWidth',1) 

        grid on; hold on; xlabel('f [Hz]'); ylabel('dB'); 

        if ii==3 

            ylabel('20log10(A)') 

        else 

            xlim([0,25e3]); ylim([0,80]) 

        end   

        %         title([s.Channels(ii).Device.Model(1:3) '...' ' Ch ' s.Channels(ii).ID... 

        %             ]); 

        title(['FFT ' title_str(ii)]); 

        x1=frequency(index(1,ii)); y1=amp_desc(1,ii); 

        x2=frequency(index(2,ii)); y2=amp_desc(2,ii); 

        str1 = ['(',num2str(x1),'Hz,',num2str(y1),')']; 

        str2=[' THD=',num2str(THD_percentage(ii)),'%']; 

        p1 = plot(x1,y1,'o','MarkerFaceColor','red'); 

        %p2= plot(x2,y2,'o','MarkerFaceColor','red'); 

        text(x1+10,y1,str1,'FontSize',8,'VerticalAlignment','top','Margin',0.1); 

        text(0.99,0.01,str2,'FontSize',10,'VerticalAlignment','bottom',... 

            'HorizontalAlignment','right','Margin',0.1,'Units','normalized',... 

            'BackgroundColor','w'); 

    end 

    hold off 

end 

%% 

%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%% plot thd 

THD_percentage=zeros(3,40); 

for nn=1:40 

    outputDurationInSeconds = test_duration+0.1; 

    test_time = (1/Fs_out)*[1:round(outputDurationInSeconds*Fs_out)]'; 

    % % Generate sinewave test tone. 

    fc = 500*nn; A = 3; 

    test_signal = A*sin(2*pi*fc*test_time); 

    %%% outputting test_signal in background and getting input in foreground 

    queueOutputData(s2, repmat(test_signal, 1, number_of_output_channel));% repmat to match output data to channel 

    lh=addlistener(s2,'DataRequired', ... 

        @(src,event) src.queueOutputData(s2, repmat(y, 1, number_of_output_channel)));% add listener for output 

    startBackground(s2); 

    %pause(0.1); % matlab has an internal delay,pause to avoid getting useless signal input 

    [data,t] = startForeground(s); 

    delete (lh) % delete listener 

    %%% 

    %%% bandpass filter 

    d = fdesign.bandpass('N,F3dB1,F3dB2',10,200,60e3,Fs); 

    Hd = design(d,'butter'); 

    %%%do this if it is differential input. 

    number_of_input_channel=3; 

    data_1=0.5*(data(:,1)-data(:,2)); 

    data_2=0.5*(data(:,3)-data(:,4)); 

    data=[data_1 data_2 data(:,5)]; 

    data2=data; 

    for ii=1:number_of_input_channel 

        data(:,ii)=filter(Hd,data2(:,ii)); 

    end 

    %%% 

    [ ~,~, THD_percentage(:,nn)] = fft_hz( data,Fs); pause(0.1) 

end 

figure(); 

plot(500*(1:40),thd1,'k','LineWidth',1);hold all 

plot(500*(1:40),thd2,'b','LineWidth',1) 

xlabel('[Hz]');ylabel('THD [%]');title('THD');grid on; 

legend('Position 1','Position 2'); 

%plottool 

% thd2=THD_percentage(1,:); 

% load thd1 

%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%% 

%% 

%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%% frequency sweep 

test_duration=20; 

s.DurationInSeconds = test_duration; 

outputDurationInSeconds = test_duration+0.1; 

test_time = (1/Fs_out)*[1:round(outputDurationInSeconds*Fs_out)]'; 

% Generate linear sweep 

test_signal = 3*chirp(test_time,10,0.8*test_duration,20000); 

% Start @ 15Hz,cross 20000Hz at t=0.8*test_duration [sec] 

% outputting test_signal in background and getting input in foreground 

queueOutputData(s2, repmat(test_signal, 1, number_of_output_channel));% repmat to match output data to channel 

lh=addlistener(s2,'DataRequired', ... 

    @(src,event) src.queueOutputData(s2, repmat(y, 1, number_of_output_channel)));% add listener for output 

startBackground(s2); 

%pause(0.1); % matlab has an internal delay,pause to avoid getting useless signal input 

[data,t] = startForeground(s); 

delete (lh) % delete listener 

% bandpass filter 

d = fdesign.bandpass('N,F3dB1,F3dB2',10,200,60e3,Fs); 

Hd = design(d,'butter'); 

%%%do this if it is differential input. 

number_of_input_channel=3; 

data_1=0.5*(data(:,1)-data(:,2)); 

data_2=0.5*(data(:,3)-data(:,4)); 

data=[data_1 data_2 data(:,5)]; 

data2=data; 

for ii=1:number_of_input_channel 

    data(:,ii)=filter(Hd,data2(:,ii)); 

end 

title_str={'Reference 1','Reference 2','Amplifier output'}; 

%%% 

for ii=1:number_of_input_channel 

    subplot(2,2,ii); plot(t, data(:,ii)); title(title_str(ii)); 

end 

[ frequency,amplitude, THD_percentage] = fft_hz( data,Fs); 

amplitude = sgolayfilt(amplitude,7,39);% curve fitting 

sensitivity=[186,186,0]; 

amplitude=repmat(sensitivity,length(amplitude),1)+amplitude; 

frequency_plot_length=round(length(frequency)/4); 

%%% 

figure();subplot(2,1,1); 

plot(frequency(1:frequency_plot_length),amplitude(1:frequency_plot_length,1),'k','LineWidth',1) 

hold on;grid on; 

plot(frequency(1:frequency_plot_length),amplitude(1:frequency_plot_length,2),'b','LineWidth',1) 

xlim([0,25e3]);ylim([70,120]);title('Reference microphones'); 

legend('Position 1','Position 2');xlabel('f [Hz]'); ylabel('dB'); 

subplot(2,1,2); plot(frequency(),amplitude(:,3),'k','LineWidth',1) 

title('Amplifier output');grid on;xlabel('f [Hz]');ylabel('20log10(A)') 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%