specscopepp

0001 function outdata=specscopepp(indata)
0002 
0003 global acq;
0004 
0005 h=hamming(5);
0006 mins=5e-008;
0007 maxs=1e-004;
0008 
0009 % record and plot audio spectrogram
0010 %
0011 % Usage: outdata=specscope(indata)
0012 %
0013 %  Input: indata (optional)
0014 %    Displays a recorded piece of data, if an argument is passed
0015 %    Otherwise displays audio data from an attached microphone
0016 %
0017 %  Output: outdata (optional)
0018 %    If present, will return up to 10 minutes
0019 %    of captured audio data.
0020 %
0021 %  Note: Parameters such as sampling frequency, number of tapers
0022 %    and display refresh rate may be set below if desired.
0023 %    You can also acquire data from a national instruments card.
0024 %
0025 
0026 close all;
0027 
0028 %=======================================================================
0029 % Check toolboxes
0030 %=======================================================================
0031 % Check for toolboxes
0032 if not(exist('analoginput','file'));
0033     fprintf('You need to install the DAQ toolbox first\n');
0034     return
0035 end
0036 if not(exist('mtspecgramc','file'));
0037     fprintf('You need to install the Chronux toolbox first from http://chronux.org/\n');
0038     return
0039 end
0040 
0041 %=======================================================================
0042 % Set parameters
0043 %=======================================================================
0044 
0045 
0046 % Set defaults
0047 acq.params.Fs = 44100;
0048 acq.pause = 0;
0049 acq.skips = 0;
0050 acq.stop = 0;
0051 acq.restart = 0;
0052 acq.plot_frequency = 20;
0053 acq.samples_acquired = 0;
0054 acq.spectra = [];
0055 acq.times = [];
0056 defaults
0057 audio_instr;
0058 fig=create_ui;
0059 
0060 %=======================================================================
0061 % Check arguments, start DAQ
0062 %=======================================================================
0063 
0064 if nargout 
0065    % save up to ten minutes data, preallocated...
0066     fprintf('Pre-allocating memory for data save - please be patient.\n');
0067    outdata=zeros( (acq.params.Fs * 60 * 10), 1 ); 
0068 end    
0069 
0070 if nargin == 1;
0071     acq.indata = indata;
0072     acq.live = 0;
0073 else
0074     % Create and set up and start analog input daq
0075     input=1;
0076     if input==1;
0077         acq.ai = analoginput('winsound');
0078         addchannel( acq.ai, 1 );
0079     else
0080         acq.ai = analoginput('nidaq', 1);
0081         addchannel(acq.ai, 0);
0082         set(acq.ai,'InputType','SingleEnded')
0083         set(acq.ai,'TransferMode','Interrupts')
0084         set(acq.ai,'TriggerType','Manual');
0085     end
0086     set( acq.ai, 'SampleRate', acq.params.Fs )
0087     acq.params.Fs = get( acq.ai, 'SampleRate' );
0088     acq.samples_per_frame = acq.params.Fs / acq.plot_frequency;
0089     set( acq.ai, 'SamplesPerTrigger', inf )
0090     start(acq.ai)
0091     acq.live = 1;
0092 
0093     if input==2;
0094         trigger(acq.ai);
0095     end
0096 end
0097 
0098 %=======================================================================
0099 % The scope main loop
0100 %=======================================================================
0101 
0102 acq.t0=clock;
0103 acq.tn=clock;
0104 
0105 % Loop over frames to acquire and display
0106 while 1;
0107 
0108     % Check for quit signal
0109     if acq.stop;
0110         break;
0111     end
0112     
0113     % Calculate times
0114     calctime = acq.samples_acquired / acq.params.Fs;
0115     acq.samples_acquired = acq.samples_acquired + acq.samples_per_frame;
0116     acq.t1 = clock;
0117     elapsed = etime(acq.t1,acq.t0);
0118 
0119     % Get a small snippet of data
0120     if ( acq.live )
0121        data = getdata( acq.ai, acq.samples_per_frame );
0122     else
0123       while elapsed < acq.samples_acquired / acq.params.Fs
0124         pause( acq.samples_acquired / (acq.params.Fs) - elapsed );
0125         acq.t1=clock;
0126         elapsed = etime(acq.t1,acq.t0);
0127       end
0128       if acq.samples_acquired + 2 * acq.samples_per_frame >= length( acq.indata )
0129         acq.stop=1;
0130       end
0131       data = acq.indata(floor(acq.samples_acquired+1):floor(acq.samples_per_frame+acq.samples_acquired));
0132     end
0133 
0134     if nargout 
0135         outdata(floor(acq.samples_acquired+1):floor(acq.samples_acquired+length(data))) = data(:);
0136     end
0137 
0138     if acq.restart;
0139        acq.restart = 0;
0140        acq.spectra = [];
0141        acq.times = [];
0142     end
0143 
0144     data=conv(h,data);
0145     data=sign(data-acq.threshold);
0146 
0147     % Calculate spectrogram of data snippet
0148     if acq.deriv
0149       [s, t, f] = mtspecgramc(abs(diff(data)), acq.moving_window, acq.params );
0150     else
0151       [s, t, f] = mtspecgramc(diff(data), acq.moving_window, acq.params );
0152     end
0153     
0154     % Add new spectra to that already calculated
0155     acq.times = [acq.times t+calctime];
0156     if acq.log
0157         acq.spectra = [acq.spectra log(s')];
0158     else
0159         acq.spectra = [acq.spectra s'];
0160     end                
0161 
0162     % Remove old spectra once window reaches desired size
0163     while acq.times(1,end) - acq.times(1,1) > acq.display_size;
0164 
0165         % Ring buffer!
0166         y = length(t);
0167         acq.times(:,1:y) = [];
0168         acq.spectra(:,1:y) = [];
0169     end
0170 
0171     % Only plot if display is keeping up with real time and not paused
0172     show_plot=1;
0173     if nargin==0
0174        if get(acq.ai, 'SamplesAvailable' ) > 10 * acq.samples_per_frame && acq.pause==0
0175       show_plot=0;
0176        end
0177     else
0178       if elapsed > calctime + 0.5
0179     show_plot=0;
0180       end
0181     end
0182 
0183     if acq.pause
0184        show_plot=0;
0185     end
0186     if show_plot
0187         
0188         if  acq.bgsub
0189             acq.mean_spectra = mean( acq.spectra, 2 );
0190         end
0191         
0192         % Normalize until full screen passes by if requested
0193         if acq.normalize>=1;
0194             if acq.normalize==1
0195                 acq.tn=clock;
0196                 acq.normalize=2;
0197             end
0198             if etime(clock,acq.tn)>1.25*acq.display_size
0199                 acq.normalize=0;
0200             end
0201             mins = min(min(acq.spectra));
0202             maxs = max(max(acq.spectra));
0203         end
0204 
0205         % Scale the spectra based upon current offset and scale
0206         if acq.bgsub
0207            scaled_spectra = acq.offset + ( acq.scale ) * ( acq.spectra - repmat( acq.mean_spectra, [1,size(acq.spectra,2)]) ) / ( maxs - mins ); 
0208         else
0209             scaled_spectra = acq.offset + acq.scale * ( acq.spectra - mins ) / ( maxs - mins );      
0210         end
0211 
0212     
0213         % Draw the image to the display
0214         image( acq.times, f, scaled_spectra, 'parent', acq.ax1 ); axis([acq.ax1],'xy');
0215         drawnow;
0216  
0217     else
0218         % Keep track of skipped displays
0219         acq.skips = acq.skips + 1;
0220     end
0221 
0222 end
0223 
0224 %=======================================================================
0225 % Clean up
0226 %=======================================================================
0227 
0228 acq.t1=clock;
0229 elapsed = etime(acq.t1,acq.t0);
0230 fprintf( 'Elapsed time %f seconds\n', elapsed );
0231 
0232 % Warn if many skips were encountered
0233 if acq.skips > 5;
0234     fprintf( '\nWARNING:\nThis program skipped plotting %d times to keep pace.\n', acq.skips )
0235     fprintf( 'Run again without keyboard interaction or changing the figure size.\n' )
0236     fprintf( 'If this message reappears you should reduce the plot frequency parameter.\n\n' );
0237 end
0238 
0239 % Clean up the analoginput object
0240 if acq.live
0241   stop(acq.ai);delete( acq.ai );clear acq.ai;
0242 end
0243 
0244 % Clean up the figure
0245 delete(fig);
0246 delete(gcf);
0247 
0248 if nargout 
0249    % save up to ten minutes data, preallocated...
0250     fprintf('Saving output data\n');
0251    outdata=outdata(1:floor(acq.samples_acquired));
0252 end  
0253 
0254 return;
0255 
0256 %
0257 %
0258 %=======================================================================
0259 % Functions called
0260 %=======================================================================
0261 %
0262 %
0263 
0264 
0265 %=======================================================================
0266 % Handle Keypresses
0267 %=======================================================================
0268 
0269 % Handle figure window keypress events
0270 function keypress(varargin)
0271 
0272 global acq;
0273 keypressed=get(gcf,'CurrentCharacter');
0274 
0275 % ignore raw control, shift, alt keys
0276 if keypressed;
0277     
0278     % Offset changes
0279     increment=1;
0280     if strcmp( keypressed, 'l');
0281     acq.offset = acq.offset - increment;
0282     elseif strcmp( keypressed, 'o');
0283         acq.offset = acq.offset + increment;
0284 
0285     % Scale changes
0286     elseif strcmp( keypressed, 'x');
0287         acq.scale = acq.scale - increment;
0288     elseif strcmp( keypressed, 's');
0289         acq.scale = acq.scale + increment;
0290 
0291     % Reset defaults
0292     elseif strcmp( keypressed, 'd');
0293         defaults
0294     acq.restart=1;
0295     % Normalize spectra
0296     elseif strcmp( keypressed, 'n');
0297        request_normalize
0298         
0299     % Quit
0300     elseif strcmp( keypressed, 'q');
0301     request_quit
0302 
0303     % Pause
0304     elseif strcmp( keypressed, 'p');
0305        request_pause
0306   
0307     % Help
0308     elseif strcmp( keypressed, 'h');
0309         audio_instr
0310 
0311 
0312     elseif strcmp( keypressed, 't');
0313        acq.threshold = acq.threshold + 0.01;
0314 
0315     elseif strcmp( keypressed, 'g');
0316        acq.threshold = acq.threshold - 0.01;
0317 
0318 
0319   
0320     % Change colormaps for 0-9,a-c
0321     elseif strcmp( keypressed, '0' );
0322         colormap( 'jet' );
0323     elseif strcmp( keypressed, '1' );
0324         colormap( 'bone' );
0325     elseif strcmp( keypressed, '2' );
0326         colormap( 'colorcube' );
0327     elseif strcmp( keypressed, '3' );
0328         colormap( 'cool' );
0329     elseif strcmp( keypressed, '4' );
0330         colormap( 'copper' );
0331     elseif strcmp( keypressed, '5' );
0332         colormap( 'gray' );
0333     elseif strcmp( keypressed, '6' );
0334         colormap( 'hot' );
0335     elseif strcmp( keypressed, '7' );
0336         colormap( 'hsv' );
0337     elseif strcmp( keypressed, '8' );
0338         colormap( 'autumn' );
0339     elseif strcmp( keypressed, '9' );
0340         colormap( 'pink' );
0341     elseif strcmp( keypressed, 'a' );
0342         colormap( 'spring' );
0343     elseif strcmp( keypressed, 'b' );
0344         colormap( 'summer' );
0345     elseif strcmp( keypressed, 'c' );
0346         colormap( 'winter' );        
0347     end
0348 
0349     update_display
0350     
0351 end
0352 return
0353 
0354 %=======================================================================
0355 % Defaults
0356 %=======================================================================
0357 
0358 % Reset defaults
0359 function defaults()
0360   global acq;
0361   acq.params.raw_tapers = [3 5];
0362   acq.moving_window = [0.02 0.01];
0363   acq.params.tapers=dpsschk(acq.params.raw_tapers,round(acq.params.Fs*acq.moving_window(1)),acq.params.Fs);
0364   acq.offset = 0;
0365   acq.scale = 500;
0366   acq.display_size = 3;
0367   acq.params.fpass = [50 20000];
0368   acq.deriv=1;
0369   acq.log=0;
0370   acq.bgsub = 1;
0371   acq.params.pad= 0;
0372   acq.normalize = 0;
0373   acq.threshold=0;
0374 
0375 return
0376 
0377 function update_display()
0378 global acq;
0379     set(acq.tapers_ui,'String',sprintf( '%.0f %.0f', acq.params.raw_tapers(1), acq.params.raw_tapers(2) ));
0380     set(acq.window_ui,'String',sprintf( '%.2f %.2f', acq.moving_window(1), acq.moving_window(2) ));
0381     set(acq.offset_ui,'String',sprintf( '%d', acq.offset ));
0382     set(acq.scale_ui,'String',sprintf( '%d', acq.scale ));
0383     set(acq.display_size_ui,'String',sprintf( '%.1f', acq.display_size ));
0384     set(acq.frequency_ui,'String',sprintf( '%.1f %.1f', acq.params.fpass(1), acq.params.fpass(2)  ))
0385     set(acq.derivative_ui,'Value',acq.deriv);
0386     set(acq.log_ui,'Value',acq.log);
0387     set(acq.bgsub_ui,'Value',acq.bgsub);
0388     set(acq.threshold_ui,'String',sprintf( '%.2f', acq.threshold ));
0389     return
0390 
0391 
0392 %=======================================================================
0393 % Update ui controls
0394 %=======================================================================
0395 
0396 function request_quit(varargin)
0397      global acq;
0398      acq.stop=1;
0399 return
0400 
0401 function request_pause(varargin)
0402      global acq;
0403      acq.pause = not( acq.pause );
0404 return
0405 
0406 function request_normalize(varargin)
0407     global acq;
0408         acq.normalize = 2;
0409 return
0410 
0411 function update_defaults(varargin)
0412   global acq;
0413   defaults
0414   update_display
0415   acq.restart=1;
0416 return
0417 
0418 function update_tapers(varargin)
0419      global acq;
0420      acq.params.raw_tapers = sscanf(get( gco, 'string' ),'%f %d')';
0421      acq.params.tapers=dpsschk(acq.params.raw_tapers,round(acq.params.Fs*acq.moving_window(1)),acq.params.Fs); % check tapers
0422 return
0423 
0424 function update_window(varargin)
0425      global acq;
0426      acq.moving_window = sscanf(get( gco, 'string' ),'%f %f');
0427          acq.params.tapers=dpsschk(acq.params.raw_tapers,round(acq.params.Fs*acq.moving_window(1)),acq.params.Fs);
0428      acq.restart = 1;
0429 return
0430 
0431 function update_offset(varargin)
0432      global acq;
0433      acq.offset = sscanf(get( gco, 'string' ),'%f');
0434      return
0435 
0436 function update_scale(varargin)
0437      global acq;
0438      acq.scale = sscanf(get( gco, 'string' ),'%f');
0439      return
0440 
0441 function update_display_size(varargin)
0442      global acq;
0443      acq.display_size = sscanf(get( gco, 'string' ),'%f');
0444      return
0445 
0446 function update_fpass(varargin)
0447      global acq;
0448      acq.params.fpass = sscanf(get( gco, 'string' ),'%f %f');
0449          acq.restart = 1;
0450      return
0451 
0452 function update_deriv(varargin)
0453      global acq;
0454      acq.deriv=get( gco, 'Value' );
0455      acq.normalize=1;
0456      return
0457      
0458 function update_log(varargin)
0459      global acq;
0460      acq.log=get( gco, 'Value' );
0461      acq.normalize=1;
0462      return
0463 
0464 function update_bgsub(varargin)
0465      global acq;
0466      acq.bgsub=get( gco, 'Value' );
0467      return
0468      
0469 function update_threshold(varargin)
0470      global acq;
0471      acq.threshold = sscanf(get( gco, 'string' ),'%f');
0472      return
0473      
0474 %=======================================================================
0475 % UI display
0476 %=======================================================================
0477 
0478 function fig=create_ui()
0479     global acq;
0480 
0481     bgcolor = [.7 .7 .7];
0482 
0483     % ===Create main figure==========================
0484     fig = figure('Position',centerfig(800,600),...
0485         'NumberTitle','off',...
0486         'Name','Real-time spectrogram',...
0487         'doublebuffer','on',...
0488         'HandleVisibility','on',...
0489     'Renderer', 'openGL', ...
0490     'KeyPressFcn', @keypress, ...
0491         'Color',bgcolor);
0492 
0493     acq.ax1 = axes('position', [0.05,0.1,0.9,0.85]);
0494 
0495     colormap( 'autumn' );
0496 
0497     % ===text==========
0498     uicontrol(gcf,'Style','text',...
0499         'String', 'tapers',...
0500         'Position',[225 20 45 20],...
0501         'BackgroundColor',bgcolor);
0502     uicontrol(gcf,'Style','text',...
0503         'String', 'moving win',...
0504         'Position',[300 20 70 20],...
0505         'BackgroundColor',bgcolor);
0506     uicontrol(gcf,'Style','text',...
0507         'String', 'offset',...
0508         'Position',[375 20 30 20],...
0509         'BackgroundColor',bgcolor);
0510     uicontrol(gcf,'Style','text',...
0511         'String', 'scale',...
0512         'Position',[410 20 30 20],...
0513         'BackgroundColor',bgcolor);
0514     uicontrol(gcf,'Style','text',...
0515         'String', 't axis',...
0516         'Position',[445 20 30 20],...
0517         'BackgroundColor',bgcolor);
0518     uicontrol(gcf,'Style','text',...
0519         'String', 'f axis',...
0520         'Position',[480 20 40 20],...
0521         'BackgroundColor',bgcolor);
0522     uicontrol(gcf,'Style','text',...
0523         'String', 'abs',...
0524         'Position',[550 20 35 20],...
0525         'BackgroundColor',bgcolor);
0526     uicontrol(gcf,'Style','text',...
0527         'String', 'log',...
0528         'Position',[580 20 35 20],...
0529         'BackgroundColor',bgcolor);
0530     uicontrol(gcf,'Style','text',...
0531         'String', 'bgsub',...
0532         'Position',[610 20 35 20],...
0533         'BackgroundColor',bgcolor);
0534     uicontrol(gcf,'Style','text',...
0535         'String', 'thresh',...
0536         'Position',[645 20 35 20],...
0537         'BackgroundColor',bgcolor);
0538  
0539     % ===The quit button===============================
0540     uicontrol('Style','pushbutton',...
0541         'Position',[5 5 45 20],...
0542         'String','Quit',...
0543         'Interruptible','off',...
0544         'BusyAction','cancel',...
0545         'Callback',@request_quit);
0546 
0547     % ===The pause button===============================
0548     uicontrol('Style','pushbutton',...
0549         'Position',[55 5 45 20],...
0550         'String','Pause',...
0551         'Interruptible','off',...
0552         'BusyAction','cancel',...
0553         'Callback',@request_pause);
0554 
0555     % ===The defaults button===============================
0556     uicontrol('Style','pushbutton',...
0557         'Position',[105 5 50 20],...
0558         'String','Defaults',...
0559         'Interruptible','off',...
0560         'BusyAction','cancel',...
0561         'Callback',@update_defaults);
0562 
0563     % ===The normalize button===============================
0564     uicontrol('Style','pushbutton',...
0565         'Position',[160 5 60 20],...
0566         'String','Normalize',...
0567         'Interruptible','off',...
0568         'BusyAction','cancel',...
0569         'Callback',@request_normalize );
0570 
0571     % ===Tapers============================================
0572     acq.tapers_ui = uicontrol(gcf,'Style','edit',...
0573         'String', sprintf( '%.0f %.0f', acq.params.raw_tapers(1), acq.params.raw_tapers(2) ),...
0574         'Position',[225 5 70 20],...
0575         'CallBack', @update_tapers);
0576 
0577     % ===Window============================================
0578     acq.window_ui=uicontrol(gcf,'Style','edit',...
0579         'String', sprintf( '%.2f %.2f', acq.moving_window(1), acq.moving_window(2) ),...
0580         'Position',[300 5 70 20],...
0581         'CallBack', @update_window);
0582 
0583     % ===Offset============================================
0584     acq.offset_ui = uicontrol(gcf,'Style','edit',...
0585         'String', sprintf( '%d', acq.offset ),...
0586         'Position',[375 5 30 20],...
0587         'CallBack', @update_offset);
0588 
0589     % ===Scale============================================
0590     acq.scale_ui = uicontrol(gcf,'Style','edit',...
0591         'String', sprintf( '%d', acq.scale ),...
0592         'Position',[410 5 30 20],...
0593         'CallBack', @update_scale);
0594 
0595     % ===display size======================================
0596     acq.display_size_ui = uicontrol(gcf,'Style','edit',...
0597         'String', sprintf( '%.1f', acq.display_size ),...
0598         'Position',[445 5 30 20],...
0599         'CallBack', @update_display_size);
0600 
0601     % ===frequency axis=====================================
0602     acq.frequency_ui = uicontrol(gcf,'Style','edit',...
0603         'String', sprintf( '%.1f %.1f', acq.params.fpass(1), acq.params.fpass(2)  ),...
0604         'Position',[480 5 80 20],...
0605         'CallBack', @update_fpass);
0606 
0607     % ===derivative=====================================
0608     acq.derivative_ui = uicontrol(gcf,'Style','checkbox',...
0609         'Value',acq.deriv,...
0610         'Position',[565 5 20 20],...
0611         'CallBack', @update_deriv);
0612     
0613     % ===log=====================================
0614     acq.log_ui = uicontrol(gcf,'Style','checkbox',...
0615         'Value',acq.log,...
0616         'Position',[590 5 20 20],...
0617         'CallBack', @update_log);
0618 
0619     % ===bgsub=====================================
0620     acq.bgsub_ui = uicontrol(gcf,'Style','checkbox',...
0621         'Value',acq.bgsub,...
0622         'Position',[615 5 20 20],...
0623         'CallBack', @update_bgsub);
0624 
0625     % ===threshold======================================
0626     acq.threshold_ui = uicontrol(gcf,'Style','edit',...
0627         'String', sprintf( '%.2f', acq.threshold ),...
0628         'Position',[640 5 40 20],...
0629         'CallBack', @update_threshold);
0630 
0631 return
0632 
0633 
0634 %=======================================================================
0635 % Assorted functions
0636 %=======================================================================
0637 
0638 function pos = centerfig(width,height)
0639 % Find the screen size in pixels
0640 screen_s = get(0,'ScreenSize');
0641 pos = [screen_s(3)/2 - width/2, screen_s(4)/2 - height/2, width, height];
0642 return
0643 
0644 
0645 function audio_instr()
0646 % Show instructions
0647 
0648   fprintf('INSTRUCTIONS:\n');
0649   fprintf('Click on figure window first to activate controls.\n')
0650   fprintf('Adjust tapers, windows, scales, offsets and axes using the gui\n');
0651   fprintf('The abs checkbox toggles abs of the data\n');
0652   fprintf('The log checkbox toggles a log of the spectrum\n');
0653   fprintf('Press d or use defaults button to reset most parameters to defaults.\n')
0654   fprintf('Press n or use normalize button to normalize spectra based upon values in current display.\n')
0655   fprintf('Press 0-9,a-c to choose a colormap (default 0).\n')
0656   fprintf('Press p to pause and unpause display.\n')
0657   fprintf('Press t and g to adjust threshold, or use offset textbox on gui.\n');
0658   fprintf('Press o and l to adjust offset, or use offset textbox on gui.\n');
0659   fprintf('Press s and x to adjust scale, or use scale textbox on gui.\n');
0660   fprintf('Press h for this message.\n')
0661   fprintf('Press q to quit, or use quit button on gui.\n\n')
0662 
0663 return
0664