%% Copyright (C) 2018 Anthony Richardson %% %% This program is free software; you can redistribute it and/or modify it under %% the terms of the GNU General Public License as published by the Free Software %% Foundation; either version 3 of the License, or (at your option) any later %% version. %% %% This program is distributed in the hope that it will be useful, but WITHOUT %% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or %% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more %% details. %% %% You should have received a copy of the GNU General Public License along with %% this program; if not, see . % Interference Removal function bugtest % Change these values to change the portion of the time waveform that % is displayed. toff is the start time. tdur is the duration % These should be GUI controls, but they don't fit in the current % design toff = 100e-3; tdur = 20e-3; % Interference level Ai = 0.3; % Interference and notch frequencies fi = 1500; fn = 1500; % Decay constant for impulse response alpha = 200; % Save initial simulation state simState('Ai', Ai); simState('fi', fi); simState('fn', fn); simState('alpha', alpha); simState('toff', toff); simState('tdur', tdur); guiState('xup', true); guiState('yup', true); guiState('zup', true); guiState('hup', true); guiState('xplay' , false); guiState('yplay' , false); guiState('zplay' , false); clf; createPanel(); createPlot(); end function createPanel pos = [.67 0.1 .30 .50 ]; p = uipanel ('position', pos, 'visible', 'off'); guiState('cpanel', p); %% BUG workaround - may be able to remove this in a %% future version of octave if (exist('OCTAVE_VERSION', 'builtin')) % The panel position seems to get messed up when the % window is resized in Octave. This was an attempt to fix it. %set(gcf, 'sizechangedfcn', @resetControlPanel); % That didn't work but this seems too. dellistener(gcf, 'position'); addlistener (gcf, 'position', {@resetControlPanel, p, pos}); end bg = uibuttongroup(p, 'position', [.0 .88 1 .12]); bx = uicontrol(bg, 'style', 'pushbutton', ... 'units', 'normalized', ... 'string', 'x', ... 'callback', @guiCallback, ... 'position', [0.45 0.1 .1 .8]); by = uicontrol(bg, 'style', 'pushbutton', ... 'units', 'normalized', ... 'string', 'y', ... 'callback', @guiCallback, ... 'position', [0.65 0.1 .1 .8]); bz = uicontrol(bg, 'style', 'pushbutton', ... 'units', 'normalized', ... 'string', 'z', ... 'callback', @guiCallback, ... 'position', [0.85 0.1 .1 .8]); uicontrol(bg, 'style', 'text', ... 'units', 'normalized', ... 'string', 'Play:', ... 'position', [0.05 0.1 .2 .8]); guiState('bx', bx); guiState('by', by); guiState('bz', bz); p = uipanel (p, 'position', [0 0 1 .88 ]); h = sliderControl(p, 'Ai', 0, ... 1, simState('Ai'), [.05 .2], ... [0.1 .9 1 .1], @guiCallback ); guiState('Aih', h); h = sliderControl(p, 'fi', 600, ... 2000, simState('fi'), [10 50], ... [0.1 .7 1 .1], @guiCallback ); guiState('fih', h); h = sliderControl(p, 'fn', 600, ... 2000, simState('fn'), [10 50], ... [0.1 .5 1 .1], @guiCallback ); guiState('fnh', h); h = sliderControl(p, 'alpha', 10, ... 500, simState('alpha'), [10 50], ... [0.1 .3 1 .1], @guiCallback ); guiState('ah', h); h = uicontrol(p, 'style', 'text', ... 'units', 'normalized', ... 'string', 'ready', ... 'backgroundcolor', [.2 .2 .2], ... 'foregroundcolor', [1 1 1], ... 'position', [0.35 0.02 .3 .1]); guiState('status', h); end % Reset GUI panel after window size change function resetControlPanel(s, e, p, pos) % Just setting the panel position does not seem to work but % maximizing it and then resetting the position does. set(p, 'visible', 'off'); set(p, 'position', [0 0 1 1]); set(p, 'position', pos); set(p, 'visible', 'on'); end function guiCallback(s, ~) g = guiState(); if (s == g.bx) guiState('xplay', true); elseif (s == g.by) guiState('yplay', true); elseif (s == g.bz) guiState('zplay', true); elseif (s == g.Aih) % Ai slider (inter amplitude) v = round(get(s, 'value')/0.05)*0.05; set(s, 'value', v); set(get(s, 'userdata'), 'string', num2str(v)); simState('Ai', v); guiState('yup', true); guiState('zup', true); elseif (s == g.fih) % fi slider (interference freq) v = round(get(s, 'value')/10)*10; set(s, 'value', v); set(get(s, 'userdata'), 'string', num2str(v)); simState('fi', v); guiState('yup', true); guiState('zup', true); elseif (s == g.fnh) % fn slider (notch freq) v = round(get(s, 'value')/10)*10; set(s, 'value', v); set(get(s, 'userdata'), 'string', num2str(v)); simState('fn', v); guiState('hup', true); guiState('zup', true); elseif (s == g.ah) % alpha slider (decay) v = round(get(s, 'value')/10)*10; set(s, 'value', v); set(get(s, 'userdata'), 'string', num2str(v)); simState('alpha', v); guiState('hup', true); guiState('zup', true); end createPlot(); end function createPlot() %%{ fs = 44100; L = 32768; % Keep x, y, and z persistent so that they can be played without % having to be recalculated persistent x if (isempty(x)) x = sin(2*pi*1000*(0:(L-1))/fs); end persistent y if(isempty(y)) y = zeros(size(x)); end persistent z if (isempty(z)) z = zeros(size(x)); end % These are persistent for performance reasons. % They only need to be computed once and never change. %persistent L %if (isempty(L)) %L = length(x) %end persistent t if (isempty(t)) t = (0:L-1)/fs; end persistent fi if (isempty(fi)) fi =(1:4000); end persistent F if (isempty(F)) F =(fi-1)*fs/L; end % These are needed to compute z and FZ persistent FH FY %} s = simState(); g = guiState(); %%{ set(g.status, 'string', 'wait'); i_sta = round(s.toff*fs)+1; i_end = i_sta + round(s.tdur*fs); nplot = i_sta:i_end; tplot = (nplot-1)*1000/fs; if(g.xup) guiState('xup', false); subplot(3,3,1); plot(tplot, x(nplot)); grid on; xlim([s.toff s.toff+s.tdur]*1000); xlabel('time (ms)'); ylabel('Amplitude'); title('Original (x)'); subplot(3,3,2); FX = fft(x); FXI = 2*abs(FX(fi))/L; semilogy(F/1000, FXI); xlim([0 5]); ylim([1e-5 1]); grid on; xlabel('freq (kHz)'); ylabel('Amplitude'); title('Original (x)'); end if(g.hup) guiState('hup', false); alpha2=s.alpha*s.alpha; wo=2*pi*s.fn; A=sqrt(4*alpha2+alpha2*alpha2/wo/wo); theta=atan(s.alpha/wo/2); h=-A*exp(-s.alpha*t).*cos(wo*t+theta); FH=fft(h)/fs; subplot(3,3,3); plot(t(1:round(L/15))*1000, h(1:round(L/15))); grid on; %xlim([0 tplotmax]); xlabel('time (ms)'); ylabel('Amplitude'); title('Impulse Response'); end if(g.yup) % 'redraw' guiState('yup', false); y=x+s.Ai *cos(2*pi*s.fi*t); subplot(3,3,4); plot(tplot, y(nplot)); grid on; xlim([s.toff s.toff+s.tdur]*1000); xlabel('time (ms)'); ylabel('Amplitude'); title('Orig + Inter (y)'); subplot(3, 3, 5) FY = fft(y); FYI = 2*abs(FY(fi))/L; semilogy(F/1000, FYI); xlim([0 5]); ylim([1e-5 1]); grid on; xlabel('freq (kHz)'); ylabel('Amplitude'); title('Orig + Inter (y)'); end if(g.zup) guiState('zup', false); subplot(3, 3, 7) z=real(ifft(FY.*FH)); z=y+z; plot(tplot, z(nplot)); grid on; xlim([s.toff s.toff+s.tdur]*1000); xlabel('time (ms)'); ylabel('Amplitude'); title('Filtered (z)'); subplot(3, 3, 8) FZ = fft(z); FZI = 2*abs(FZ(fi))/L; FHI=abs(1+FH(fi)); semilogy(F/1000, FZI); hold on; semilogy(F/1000, FHI, 'r'); hold off; xlim([0 5]); ylim([1e-5 1]); grid on; xlabel('freq (kHz)'); ylabel('Amplitude'); title('Filtered (z)'); end if(g.xplay) xplayer = audioplayer(x, fs); playblocking(xplayer); guiState('xplay', false); end if(g.yplay) yplayer = audioplayer(y, fs); playblocking(yplayer); guiState('yplay', false); end if(g.zplay) zplayer = audioplayer(z, fs); playblocking(zplayer); guiState('zplay', false); end %} set(g.cpanel, 'visible', 'on'); set(g.status, 'string', 'ready'); end function v = guiState(key, val) persistent state; if (nargin == 0) v = state; elseif (nargin == 1) v = state.(key); elseif (nargin == 2) state.(key) = val; v = val; end end function v = simState(key, val) persistent state; if (nargin == 0) % Return the entire struct v = state; elseif (nargin == 1) % Return a particular field from the struct v = state.(key); elseif (nargin == 2) % Set a key value (key must be a string) state.(key) = val; v = val; end end % Create a slider with decorations (label and min, max, and current values) function h = sliderControl(parent, label, min, max, val, step, pos, cback) % sliderstep = [minstep maxstep] % where minstep is fraction of range on arrow click % and maxstep is fraction of range on trough click h = uicontrol(parent, 'style' ,'slider',... 'units', 'normalized',... 'position', [pos(1)+0.1*pos(3) pos(2) 0.7*pos(3) pos(4)],... 'min', min, 'max', max,... 'sliderstep', [step(1)/(max-min) step(2)/(max-min)],... 'value', val, 'callback', cback); uicontrol(parent, 'style', 'text',... 'units', 'normalized',... 'position', [pos(1) pos(2) 0.1*pos(3) pos(4)],... 'String', num2str(min)); %addproperty('my_slider_min', h, 'handle', hc); uicontrol(parent, 'style', 'text',... 'units', 'normalized',... 'position', [pos(1)+0.8*pos(3) pos(2) 0.1*pos(3) pos(4)],... 'String', num2str(max)); %addproperty('my_slider_max', h, 'handle', hc); uicontrol(parent, 'style', 'text',... 'units', 'normalized',... 'position', [pos(1)+0.15*pos(3) pos(2)-pos(4) 0.4*pos(3) pos(4)],... 'horizontalalignment','left',... 'String', label); hc = uicontrol(parent, 'style', 'text',... 'units', 'normalized',... 'position', [pos(1)+.5*pos(3) pos(2)-pos(4) 0.3*pos(3) pos(4)],... 'String', num2str(val)); %addproperty('my_slider_val', h, 'handle', hc); set(h, 'userdata', hc); end