c09cc_demo

NAG Toolbox for MATLAB demo showing signal noise reduction using wavelets
Description
The GUI

function c09cc_demo

NAG Toolbox for MATLAB demo showing signal noise reduction using wavelets

Description

This is to demostrate the use of:

the 1D multi-level discrete wavelet transform routine c09cc (to calculate the sets of wavelet coefficients at different levels) and
its corresponding inverse routine c09cb (which transform back to give filtered data set).

Filtering here refers to the manipulation of the wavelet coefficients.
In terms of noise reduction this usually means reducing the magnitude
of the coefficients at the lowest levels.

The GUI

Various control values must be selected in turn. The stages are:

A one-dimensional signal, in the range [0,2], is plotted in blue.
A "Noise" slider appears. Moving the slider adds random noise of varying amplitude to the original signal. The signal + noise is replotted in orange.
A slider appears on the right hand side: "Lowest Level Filter"; This controls the factor to multiply wavelet coefficients at the lowest level by. This represents highest frequencies in the signal, so a value close to zero will remove most of the highest frequency oscillations.
When a lowest level filter factor is selected, another slider appears: "Second Lowest Level". Choosing a small value for this (close to zero) removes most of the second highest range of oscillations. A value of 1 will leave these coefficients as they are (unfiltered).
The final right hand side slider appears (for the third lowest level of wavelet coefficients). A value must be selected for this also.
The Wavelet transform is now chosen using two list menus on the left hand side: "End extension" selection and "Wavelet" selection.
When an end-extension and wavelet are selected the "Filter" button produces a third curve which is the result of applying the inverse wavelet transform to filtered wavelet coefficients. For three filter levels close to zero you should see a reduction in the noise in the signal and a gold curve approximating the original blue data curve.

  Different wavelets will perform differently .

% create structure of handles
global handles;

%            Initial Plot Properties

% Set parameters for size and position of plot.
scrsz = get(0, 'ScreenSize');
figure_handle = figure('Name', 'NAG routines c09cc and c09cd - Multi-level Wavelet Transform and Inverse', ...
		       'Position', [scrsz(1)+100 scrsz(2)+20 scrsz(3)*0.8 scrsz(4)*0.8], ...
		       'NumberTitle', 'off');
% Use double buffering to avoid flashing of animation
set(figure_handle, 'Renderer', 'Painters', 'DoubleBuffer', 'On');
axmax = 2.0;
axmin = 0.0;
aymax = 1.5;
aymin = -0.5;
axis([axmin axmax aymin aymax]);

% Add some data to the Gui properties via handles. %

handles = guihandles(figure_handle);
handles.axmax = axmax;
handles.axmin = axmin;
handles.aymax = aymax;
handles.aymin = aymin;
handles.fontSize = 11.5;
handles.axes = gca;   % Get current axes - so we can reuse in each plot() command.

handles.n = nag_int(200);
handles.t = rand(handles.n,1);
handles.x = handles.t;
handles.y = handles.t;
handles.filt = handles.t;
handles.noise = handles.t;
handles.amp=0.2;
handles.wavnam='HAAR';
handles.mode='Z';
handles.wtrans = 'Multilevel';

% Initial un-noised data
n = double(handles.n);
for k = 1:n %#ok<FXUP>
  xx = (2*k-1)/n;
  handles.t(k) = xx;
  handles.x(k) = 0.5 - abs(1-xx) + abs(0.5 - 2*xx + floor(2.0*xx));
end

handles.z = handles.x;

% Initialize the generator to a non-repeatable sequence
% genid and subid identify the base generator
genid = nag_int(1);
subid =  nag_int(1);
[handles.state, ~] = g05kg(genid, subid);
guidata(figure_handle,handles);

% plot the un-noised curve

set(handles.axes, 'FontSize', handles.fontSize);
xlabel(handles.axes,'x');
title(handles.axes,'Using NAG routines c09ca and c09cb for wavelet filtering', ...
      'FontSize', 14);
ylabel(handles.axes, 'Pure Signal', 'FontSize', 14);
hold(handles.axes);
plot(handles.axes, handles.t, handles.z);
hold(handles.axes);

guidata(figure_handle,handles);


% Control Buttions

% Control Button selecting noise amplitude and plotting noised data

h_amp_text=uicontrol(figure_handle,'Style', 'text', 'String', 'Noise Amplititude',...
          'HorizontalAlignment', 'left', 'Position', [20 330 100 15]);
h_slider=uicontrol(figure_handle,'Style', 'slider', 'max', 1.0, ...
          'Position', [20 300 100 30], 'TooltipString', ...
          'Select maximum amplitude of background noise to add', ...
	  'Callback', @setamp);
handles.amp_text=h_amp_text;
handles.slider=h_slider;
%             Control Button for Wavelet Type
guidata(figure_handle,handles);
h_wave_text=uicontrol(figure_handle,'Style', 'text', 'String', 'Wavelet', ...
	  'HorizontalAlignment', 'left', 'Visible','off',...
          'Position', [20 130 100 15]);
h_wave=uicontrol(figure_handle,'Style', 'popup', 'String', ...
	  'Haar|DB2|DB3|DB4|DB5|DB6|DB7|DB8|DB9|DB10|BIOR1.1|BIOR1.3|BIOR1.5|BIOR2.2|BIOR2.4|BIOR2.6|BIOR2.8|BIOR3.1|BIOR3.3|BIOR3.5|BIOR3.7', ...
          'Position', [20 100 100 30], 'Visible','off','TooltipString', 'Select Wavelet');
handles.wave_text=h_wave_text;
handles.wave=h_wave;
guidata(figure_handle,handles);
%            Control Button for End Extension Method

h_end_text=uicontrol(figure_handle,'Style', 'text', 'String', 'End Extension method', ...
	  'HorizontalAlignment', 'left', 'Visible','off',...
          'Position', [20 230 100 15]);
h_end=uicontrol(figure_handle,'Style', 'popup', 'String',...
		'Half-point symmetric|Whole-point symmetric|Zero', ...
		'Position', [20 200 100 30],  'Visible','off','TooltipString', 'Select mode');
handles.end_text=h_end_text;
handles.end=h_end;
guidata(figure_handle,handles);
% Control Buttons selecting filter levels

h_f1_t=uicontrol(figure_handle,'Style', 'text', 'String', 'Lowest Level Filter',...
          'HorizontalAlignment', 'left', 'Units','normalized',...
	  'Position', [0.91 0.83 0.04 0.05],'Visible','off');
h_f1_text=uicontrol(figure_handle,'Style', 'edit', 'Units','normalized',...
			  'Position', [0.91 0.79 0.05 0.03], 'TooltipString',...
			  'you can enter a number between 0.0 and 1.0 here',  'Visible','off',...
			  'Callback',@slider_editText);
h_filt1=uicontrol(figure_handle,'Style', 'slider', 'max', 1.0, ...
          'Units','normalized','Position', [0.97 0.78 0.02 0.1], 'TooltipString', ...
          'Factor to apply to lowest level coefficients','Visible','off',...
          'Callback',@slider_filt1);
handles.f1_t=h_f1_t;
handles.f1_text=h_f1_text;
handles.filt1=h_filt1;
guidata(figure_handle,handles);
h_f2_t=uicontrol(figure_handle,'Style', 'text', 'String', 'Second Lowest Level Filter',...
          'HorizontalAlignment', 'left', 'Units','normalized',...
	  'Position', [0.91 0.73 0.04 0.06],'Visible','off');
h_f2_text=uicontrol(figure_handle,'Style', 'edit', 'Units','normalized',...
			  'Position', [0.91 0.69 0.05 0.03], 'TooltipString',...
			  'you can enter a number between 0.0 and 1.0 here',  'Visible','off',...
			  'Callback',@slider_editText2);
h_filt2=uicontrol(figure_handle,'Style', 'slider', 'max', 1.0, ...
          'Units','normalized','Position', [0.97 0.68 0.02 0.1], 'TooltipString', ...
          'Factor to apply to second lowest level coefficients','Visible','off',...
          'Callback',@slider_filt2);
handles.f2_t=h_f2_t;
handles.f2_text=h_f2_text;
handles.filt2=h_filt2;
guidata(figure_handle,handles);
h_f3_t=uicontrol(figure_handle,'Style', 'text', 'String', 'Third Lowest Level Filter',...
          'HorizontalAlignment', 'left', 'Units','normalized',...
	  'Position', [0.91 0.62 0.04 0.06],'Visible','off');
h_f3_text=uicontrol(figure_handle,'Style', 'edit', 'Units','normalized',...
			  'Position', [0.91 0.58 0.05 0.03], 'TooltipString',...
			  'you can enter a number between 0.0 and 1.0 here',  'Visible','off',...
			  'Callback',@slider_editText3);
h_filt3=uicontrol(figure_handle,'Style', 'slider', 'max', 1.0, ...
          'Units','normalized','Position', [0.97 0.57 0.02 0.1], 'TooltipString', ...
          'Factor to apply to third lowest level coefficients','Visible','off',...
          'Callback',@slider_filt3);
handles.f3_t=h_f3_t;
handles.f3_text=h_f3_text;
handles.filt3=h_filt3;
guidata(figure_handle,handles);
% Control Button to Apply Wavelet Transforms to filter lower levels

h_filter=uicontrol(figure_handle,'Position', [20 20 40 40], 'String', 'Filter',...
                    'Visible','off', 'Callback', @apply_filter);
handles.filter=h_filter;
% save the structure
guidata(figure_handle,handles);

function setamp(hObject,eventData)
  % Get Gui data
  myhandles = guidata(gcbo);
  % Get Amplitude from slider
  myhandles.amp = get(hObject,'Value');
  set(h_f1_t,'Visible','on');
  set(h_f1_text,'Visible','on');
  set(h_filt1,'Visible','on');
  myhandles.f1_t=h_f1_t;
  myhandles.f1_text=h_f1_text;
  myhandles.filt1=h_filt1;
  % generate some random variates to add noise
  [myhandles.state, myhandles.noise, ~] = g05sa(myhandles.n, myhandles.state);

  % Clear Axes;
  cla(myhandles.axes);

  % Add Noise to Profile in z
  for j = 1:myhandles.n
    myhandles.x(j) = myhandles.amp*(myhandles.noise(j)-0.5) +  myhandles.z(j);
  end

  % Plot Noised data on top of original data.
  xlabel(myhandles.axes, 'x');
  ylabel(myhandles.axes, 'Noisy signal');
  tt = sprintf('Noisy signal with noise amplitude = %5.3f', myhandles.amp);
  title(myhandles.axes, tt, 'FontSize', 14);
  axis(myhandles.axes,[myhandles.axmin myhandles.axmax myhandles.aymin myhandles.aymax]);
  hold(myhandles.axes);
  plot(myhandles.axes,myhandles.t,myhandles.z,myhandles.t,myhandles.x);
  hold(myhandles.axes);

  % Save amplitude in Gui data.
  guidata(gcbo,myhandles);
end


function apply_filter(hObject, eventData)
  % Get GUI data.
  myhandles = guidata(gcbo);
  noise = myhandles.noise;
  n = myhandles.n;
  z = myhandles.z;
  t = myhandles.t;


  % Set up private data
  y = rand(n,1);
  x = y;
  filt = y;
  amp = get(h_slider,'Value');
  wtrans = 'Multilevel';
  val = get(h_end,'Value');
  if val == 1
    mode = 'H';
  elseif val == 2
    mode = 'W';
  elseif val == 3
    mode = 'Z';
  end
  val = get(h_wave,'Value');
  if val == 1
    wavnam = 'HAAR';
  elseif val == 2
    wavnam = 'DB2';
  elseif val == 3
    wavnam = 'DB3';
  elseif val == 4
    wavnam = 'DB4';
  elseif val == 5
    wavnam = 'DB5';
  elseif val == 6
    wavnam = 'DB6';
  elseif val == 7
    wavnam = 'DB7';
  elseif val == 8
    wavnam = 'DB8';
  elseif val == 9
    wavnam = 'DB9';
  elseif val == 10
    wavnam = 'DB10';
  elseif val == 11
    wavnam = 'BIOR1.1';
  elseif val == 12
    wavnam = 'BIOR1.3';
  elseif val == 13
    wavnam = 'BIOR1.5';
  elseif val == 14
    wavnam = 'BIOR2.2';
  elseif val == 15
    wavnam = 'BIOR2.4';
  elseif val == 16
    wavnam = 'BIOR2.6';
  elseif val == 17
    wavnam = 'BIOR2.8';
  elseif val == 18
    wavnam = 'BIOR3.1';
  elseif val == 19
    wavnam = 'BIOR3.3';
  elseif val == 20
    wavnam = 'BIOR3.5';
  elseif val == 21
    wavnam = 'BIOR3.7';
  end

  % Set up wavelet and filter

  % Query wavelet filter dimensions

  [nwl, ~, nwc, icomm, ~] = ...
     c09aa(wavnam, wtrans, mode, n);


  filt(nwl+1) = get(h_filt1,'Value');
  filt(nwl) = get(h_filt2,'Value');
  filt(nwl-1) = get(h_filt3,'Value');

  for ii = 1:n
    x(ii) = amp*(noise(ii)-0.5) +  z(ii);
  end

  [c, dwtlev, icomm, ~] = c09cc(x, nwc, nwl, icomm);

  nnz = 0;
  for ii = 1:nwl + 1
    nnz = nnz + dwtlev(ii);
  end

  nlvout = nag_int(3);

  ilev = 1;
  for k = 1:nwl + 1 - nlvout %#ok<*FXUP>
    ilev = ilev + dwtlev(k);
  end
  for k = nwl + 2 - nlvout:nwl + 1
    for l = ilev:ilev + dwtlev(k) - 1
      c(l) = c(l)*filt(k);
    end
    ilev = ilev + dwtlev(k);
  end

  % Perform inverse transform
  [y, ~] = c09cd(nwl, c, n, icomm);

  % Clear Axes.
  cla(myhandles.axes);
  xlabel(myhandles.axes, 'x');
  ylabel(myhandles.axes, '(filtered) signal');
  tt = sprintf('Signal using: wavelet = %8s, end condition = %s, levels of filtering = %3d, noise amplitude = %5.3f', ...
	           wavnam, mode, nlvout, amp);
  title (myhandles.axes,tt, 'FontSize', 14);
  axis(myhandles.axes,[myhandles.axmin myhandles.axmax myhandles.aymin myhandles.aymax]);

  hold(myhandles.axes);
  plot(myhandles.axes, t,z,t,x,t,y);
  hold(myhandles.axes);
  guidata(gcbo,myhandles);

  % wait for the Start button to be pressed before continuing
end


function button_handler(hObject, eventData)
  global quit;
  label = get(hObject, 'String');
  if strcmp(label, 'Close') || strcmp(label, 'Quit')
    close;
    quit = 1;
  elseif strcmp(label, 'Pause')
    set(hObject, 'String', 'Resume');
    uiwait;
  elseif strcmp(label, 'Start') || strcmp(label, 'Resume')
    set(hObject, 'String', 'Pause');
    uiresume;
  elseif strcmp(label, 'Next')
    uiresume;
  end
end


function slider_editText(hObject, eventData)
  myhandles = guidata(gcbo);
  % get the string for the editText component
  sliderValue = get(hObject,'String');
  %convert from string to number if possible, otherwise returns empty
  sliderValue = str2num(sliderValue);

  % if user inputs something is not a number, or if the input is less than 0
  % or greater than 100, then the slider value defaults to 0
  if (isempty(sliderValue) || sliderValue < 0.0 || sliderValue > 1.0)
    set(h_filt1,'Value',1.0);
  else
    set(h_filt1,'Value',sliderValue);
  end
  set(h_f1_text,'String', num2str(sliderValue));
  set(h_f2_t,'Visible','on');
  set(h_f2_text,'Visible','on');
  set(h_filt2,'Visible','on');
  myhandles.filt1=h_filt1;
  myhandles.f1_text=h_f1_text;
  myhandles.f2_t=h_f2_t;
  myhandles.f2_text=h_f2_text;
  myhandles.filt2=h_filt2;
  guidata(gcbo,myhandles);
end

function slider_filt1(hObject, eventData)
  myhandles = guidata(gcbo);

  %obtains the slider value from the slider component
  sliderValue = get(hObject,'Value');
  %puts the slider value into the edit text component
  set(h_f1_text,'String', num2str(sliderValue));
  set(h_f2_t,'Visible','on');
  set(h_f2_text,'Visible','on');
  set(h_filt2,'Visible','on');
  myhandles.f1_text=h_f1_text;
  myhandles.f2_t=h_f2_t;
  myhandles.f2_text=h_f2_text;
  myhandles.filt2=h_filt2;
  guidata(gcbo,myhandles);
end


function slider_editText2(hObject, eventData)
  myhandles = guidata(gcbo);
  % get the string for the editText component
  sliderValue = get(hObject,'String');
  %convert from string to number if possible, otherwise returns empty
  sliderValue = str2num(sliderValue);

  % if user inputs something is not a number, or if the input is less than 0
  % or greater than 100, then the slider value defaults to 0
  if (isempty(sliderValue) || sliderValue < 0.0 || sliderValue > 1.0)
    set(h_filt2,'Value',1.0);
  else
    set(h_filt2,'Value',sliderValue);
  end
  set(h_f2_text,'String', num2str(sliderValue));
  set(h_f3_t,'Visible','on');
  set(h_f3_text,'Visible','on');
  set(h_filt3,'Visible','on');
  myhandles.filt2=h_filt2;
  myhandles.f2_text=h_f2_text;
  myhandles.f3_t=h_f3_t;
  myhandles.f3_text=h_f3_text;
  myhandles.filt3=h_filt3;
  guidata(gcbo,myhandles);
end

function slider_filt2(hObject, eventData)
  myhandles = guidata(gcbo);

  %obtains the slider value from the slider component
  sliderValue = get(hObject,'Value');
  %puts the slider value into the edit text component
  set(h_f2_text,'String', num2str(sliderValue));
  set(h_f3_t,'Visible','on');
  set(h_f3_text,'Visible','on');
  set(h_filt3,'Visible','on');
  myhandles.f2_text=h_f2_text;
  myhandles.f3_t=h_f3_t;
  myhandles.f3_text=h_f3_text;
  myhandles.filt3=h_filt3;
  guidata(gcbo,myhandles);
end

function slider_editText3(hObject, eventData)
  myhandles = guidata(gcbo);
  % get the string for the editText component
  sliderValue = get(hObject,'String');
  %convert from string to number if possible, otherwise returns empty
  sliderValue = str2num(sliderValue);

  % if user inputs something is not a number, or if the input is less than 0
  % or greater than 100, then the slider value defaults to 0
  if (isempty(sliderValue) || sliderValue < 0.0 || sliderValue > 1.0)
    set(h_filt3,'Value',1.0);
  else
    set(h_filt3,'Value',sliderValue);
  end
  set(h_f3_text,'String', num2str(sliderValue));
  set(h_wave_text,'Visible','on');
  set(h_wave,'Visible','on');
  set(h_end_text,'Visible','on');
  set(h_end,'Visible','on');
  set(h_filter,'Visible','on');
  myhandles.f3_text=h_f3_text;
  myhandles.wave_text=h_wave_text;
  myhandles.wave=h_wave;
  myhandles.end_text=h_end_text;
  myhandles.end=h_end;
  myhandles.filter=h_filter;
  guidata(gcbo,myhandles);
end

function slider_filt3(hObject, eventData)
  myhandles = guidata(gcbo);

  %obtains the slider value from the slider component
  sliderValue = get(hObject,'Value');
  %puts the slider value into the edit text component
  set(h_f3_text,'String', num2str(sliderValue));
  set(h_wave_text,'Visible','on');
  set(h_wave,'Visible','on');
  set(h_end_text,'Visible','on');
  set(h_end,'Visible','on');
  set(h_filter,'Visible','on');
  myhandles.f1_text=h_f1_text;
  myhandles.wave_text=h_wave_text;
  myhandles.wave=h_wave;
  myhandles.end_text=h_end_text;
  myhandles.end=h_end;
  myhandles.filter=h_filter;
  guidata(gcbo,myhandles);
end

end

Contents

NAG Toolbox for MATLAB demo showing signal noise reduction using wavelets

Description

The GUI