与傅里叶的故事

什么都不说,先放一个肌电手环的matlab

function PianoRecord
%   Derived from ScopeMath_Simple (available on the MATLAB File Exchange) August-21-2007 (Gautam.Vallabha@mathworks.com)
%   Version 1.0 March-19-2013 (Chris Rorden)
%   Version 2.0 2014-08-06 (GUO YI)
%   Version 2.1 2014-10-10 (GUO YI) for new protocol
%   Version 2.2 2014-11-03 (GUO YI) 加入了同步录音功能,采集肌肉电与音频信号
%   Version 2.3 2014-11-14 (GUO YI) 修改了几个bug,加入了同步录音功能,同步采集肌肉电、惯导、音频信号
%   Version 2.4 2014-12-19 (GUO YI) 添加功能,采集数据与分析
%   Version 2.5 2015-04-17 (GUO YI) tcp
global v


v.hzChoices = [250 500 1000 2000 4000]; %acquisition rate: maximum depends on hardware
v.chChoices = [1 2 3 4 5 6 7 8]; %number of channels: maximum depends on hardware (Leonardo supports 6, Teensy3 supports 14)
v.saveDataDefault = 0; %if 1 data will be saved to disk by default, if 0 data will be discarded


%默认1000HZ
v.hzDefaultIndex = 3; %e.g. if 2, 2nd option of HzChoices is the default
%默认8通道
v.chDefaultIndex = 8; %e.g. if 2, 2nd option of ChChoices is the default
v.serDefaultIndex = 0; % 0 for auto-select, 1 for simulated data, 2 for Arduino on port 1, 3 for Arduino on port 3, etc
%显示最后2s的数据
v.gGraphTotalTimeSec = 4.0; %e.g. if 1 then the last 1 second of data is displayed


v.AutoScaleAxis = true;
v.testSignal =false;


%Typically no need to edit lines below....
v.gSecPerScreenRefresh = 0.8; %e.g. if 0.1 then screen is updated 10 times per second
v.deviceType = 0; %0 for simulated data, 1 for Arduino


v.gOscHz = 1000; %current sampling rate  采样频率
% ads1294 4通道
v.gOscChannels = 6; %current number of channels we are recording
%频谱展示
v.showPower = true; %should we conduct a FFT on recent data?
v.showPowerChannel = 1;


v.gSecPerSample = [];


v.gGraphSamples = [];
%时间
v.xData = []; %time values for samples
% x: channels; y: samples  4
v.yData = []; %most recent data channels*samples
v.gUnsavedData = []; %data not yet saved to disk


%存储了所有EMG的数据
v.totalData = [];
%惯导的原始数据
v.mpuData = [];
%EU数据
v.euData = [];
%acc数据
v.accData = [];
%四元数
v.qData = [];


v.xUnits = 'Sec';
v.yUnits = 'Signal';


v.gSaveNumber = 1; %counts time since last save
v.gSaveEveryNRefreshes = 10; %eg. if 10, then we will save to disk every tenth screen update


v.gSampleNumber = 1; %only used to determine phase for simulated data
v.gSaveDataBaseFilename = []; %
v.kBPS = 460800; %baud rate for device


% GUI variables
v.hFigure = [];
v.hAxisRaw = [];%[0 xData(end) -1 1];
v.hAxesRaw = [];
v.hAxesMath = [];
v.hStartButton = [];
v.hTriggerCheck = [];
v.hSaveCheck = [];
v.hSerialPopup = []; %popup menu lists available serial ports
v.hChannelPopup = []; %popup menu lists available channels to record
v.hHzPopup = []; %popup menu lists available sampling rates


v.hModel = []; %模式,训练模式或者是测试模式 


v.acquiringData = false; %are we currently sampling data?
v.serialObj = []; %serial port
v.rawData = []; %raw data from device - data that still needs to be decoded into discrete samples


% set up a timer to periodically get the data
% 刷新屏幕的时间间隔
v.timerObj = timer('Period', v.gSecPerScreenRefresh, 'ExecutionMode', 'fixedSpacing', 'timerFcn', {@getDataFromDeviceType});


srate = 1000;
v.b2 = fir1(srate, 95/srate*2, 'low');
v.b = fir1(srate, 5/srate*2, 'high');
v.b1 = fir1(srate, [48/srate*2,52/srate*2], 'stop');




[v.fl,v.bl] = butter(8,95/srate*2, 'low');
[v.fh,v.bh] = butter(8,5/srate*2, 'high');
[v.fno,v.bno] = butter(8,[48/srate*2,52/srate*2], 'stop');


v.channel_useful = 1:6;


%修改的配置
v.testmodel = true;
v.tcpserver = 0;
%音频存储
v.audioOn = false;
v.audiofs = 44100; 
v.saveaudiofile = 'D:\研发\Bluetooth4.0 Workspace\Matlab Data Analysis\demo板肌肉电分析\Pianote raw data\';
v.saverawfile = 'D:\研发\Bluetooth4.0 Workspace\Matlab Data Analysis\demo板肌肉电分析\Pianote raw data\';
v.savefiltfile = 'D:\研发\Bluetooth4.0 Workspace\Matlab Data Analysis\demo板肌肉电分析\Pianote raw data\';


v.datestr = datestr(now,'yymmdd_HHMMSS');


v.recorder = audiorecorder(v.audiofs, 16, 1);


v.totalsamples = 0;
 v.myRecording =[];


makeGUI(); %set up user interface


    function [newdata] = hzfilter(olddata)
        global v;
     %   fprintf('1xxxxxx1 \n');
        tmpdata= double(olddata);
        %tmpdata = tmpdata/(pow2(23)-1)*2.4;
      %  fprintf('11111111  length data %d\n',length(olddata));
       % s1 = filter(v.b, 1, tmpdata);
       % s2 = filter(v.b2, 1, s1);
      %  s3 = filter(v.b1, 1, s2);
       % f1 = filter(v.fl,v.bl,tmpdata);
       %f2 = filter(v.fh,v.bh,f1);
       % f3 = filter(v.fno,v.bno,f2);
      %  newdata = s3;
        tmpdata = filtfilt(v.b,1,tmpdata');
        tmpdata = filtfilt(v.b2,1,tmpdata);
        tmpdata = filtfilt(v.b1,1,tmpdata);
      %  fprintf('11dddd111 \n');
       % tmpdata = tmpdata(:,v.channel_useful)';
     %   fprintf('11eeee\n');
        newdata = tmpdata';
    %    fprintf('122222211\n');
    %%end
    
    
    %%----从设备获取数据, 每次刷新间隔都会调用
    function getDataFromDeviceType(hObject, eventdata)
        global v;
        [newData] = getDataFromDeviceType2();
%         switch v.deviceType
%             case 1
%                 [newData] = getDataFromDeviceType1();
%             otherwise
%                 [newData] = getDataFromDeviceType0();
%         end; %switch deviceType
        if length(newData) < 1
            return
        end;
        %newData(1,:) = hzfilter(newData(1,:));
        %fprintf('-------\n');
        v.totalData =[v.totalData newData(1:6,:)];
        
       % newData = hzfilter(newData);
        v.yData = [v.yData newData]; %append new data ,新的采样 4*simples
        %fprintf('--1-----\n');
        %取新加进来的一部分数据
        v.yData = v.yData(:,length(v.yData)-v.gGraphSamples+1:end);
       % size(v.yData)
        %v.yData= hzfilter(v.yData);
       % size(v.yData)
       % fprintf('--2----\n');
        if ~isempty(v.gSaveDataBaseFilename)
            v.gUnsavedData  = [v.gUnsavedData newData]; %append new data
            v.gSaveNumber = v.gSaveNumber + 1;
            if (v.gSaveNumber > v.gSaveEveryNRefreshes)  %flush unsaved data to disk
                flushSaveDataSub();   %每隔一段时间写一次
            end; %if SaveNumber
        end %if SaveData
        % check the user closed the window while we were waiting
        % for the device to return the waveform data
        %fprintf('--3-----\n');
        if ishandle(v.hFigure),
            axes(v.hAxesRaw1);
            plot(v.xData,v.yData(1,:));
            axes(v.hAxesRaw2);
            plot(v.xData,v.yData(2,:));
            axes(v.hAxesRaw3);
            plot(v.xData,v.yData(3,:));
            axes(v.hAxesRaw4);
            plot(v.xData,v.yData(4,:));
            axes(v.hAxesRaw5);
            plot(v.xData,v.yData(5,:));
            axes(v.hAxesRaw6);
            plot(v.xData,v.yData(6,:));
            if (v.AutoScaleAxis == false)
                axis(v.hAxisRaw);
            end;
            xlabel(v.xUnits); ylabel(v.yUnits);
            if v.showPower
                axes(v.hAxesfft);
                [freq,fftdata] = powerSpectrum(length(v.yData(1,:)), v.yData(v.showPowerChannel,:));
                plot(freq, fftdata);
                xlabel('Frequency (Hz)'); ylabel('Amplitude');
            end;%showPower
        end
       % fprintf('--4-----\n');
   %end getDataFromDeviceType()


   %%-模拟数据
   function [newSamples0] = getDataFromDeviceType0()
       %create - we create precisely the same number of observations per
       %  screen refresh, so sampling rate may be approximate
       global v;
       new = round(v.gOscHz*v.gSecPerScreenRefresh);
       if (new < 1)
           new = 1;
       end
       
       for s = 1:new,
           for c=1:v.gOscChannels,f
               newSamples0(c,s) = sin(30*2*pi*(v.gSecPerSample*(s+v.gSampleNumber)) ) + randn()*0.2;
           end;
       end
       v.gSampleNumber = v.gSampleNumber + new;
       
   %end getDataFromDeviceType2()
   %%从myo设备获取数据
   function [newSamples2] = getDataFromDeviceType2()
       global v;
       packetBytes = 4 + (2 * 6); %4+2*4 = 12bytes
       count = packetBytes;
       if (v.tcpserver.BytesAvailable + length(v.rawData)  < packetBytes)
           newSamples2 = [];
           return;
       end;
       [newRawData,count] = fread(v.tcpserver,v.tcpserver.BytesAvailable,'uchar');
       %fprintf('serial bytes: %d new  \n', length(newRawData));
       v.rawData = [v.rawData; newRawData];
       if (length(v.rawData) < 1)
           return;
       end;
       [newSamples2, v.rawData] = serDecodeSub(v.rawData);
   %end getDataFromDeviceType2()
   
   %%从myo设备获取数据
   function [newSamples1] = getDataFromDeviceType1()
       global v;
       packetBytes = 4 + (2 * (v.gOscChannels-1)); %4+2*4 = 12bytes
       count = packetBytes;
       if (v.serialObj.BytesAvailable + length(v.rawData)  < packetBytes)
           newSamples1 = [];
           return;
       end;
       [newRawData,count] = fread(v.serialObj,v.serialObj.BytesAvailable,'uchar');
       %fprintf('serial bytes: %d new  \n', length(newRawData));
       v.rawData = [v.rawData; newRawData];
       if (length(v.rawData) < 1)
           return;
       end;
       [newSamples1, v.rawData] = serDecodeSub(v.rawData);
   %end getDataFromDeviceType1()




    %%---------------------------------------------------
    function [freq,fftdata] = powerSpectrum(x,y)
        n = x;
        Fs = 1000;
        freq = ((0:n-1)./n)*Fs;
        fftdata = 20*log10(abs(fft(y)));
        idx = 1:floor(length(freq)/2);
        freq = freq(idx);
        fftdata = fftdata(idx);
    %end powerSpectrum()








    %%---------------------------------------------------
    function makeGUI()
        global v;
        v.hFigure = figure('deleteFcn', {@figureCloseCallback},'name','Biomo-V3','units','pixels','position',[50, 50, 1300, 650]);
        %set(gcf,'outerposition',get(0,'screensize'));
        % if v.showPower
        %   v.hAxesRaw  = axes('position', [0.05  0.60  0.9 0.35]);
        %    title('Raw Data');
        % %   v.hAxesMath = axes('position', [0.05  0.15  0.9 0.35]);
        %   title('Processed Data');
        %  else
        v.hAxesRaw1  = axes('position', [0.05  0.89  0.7 0.09]);
       % title('channel 1  Data');
        v.hAxesRaw2  = axes('position', [0.05  0.78  0.7 0.1]);
        %title('channel 2  Data');
        v.hAxesRaw3  = axes('position', [0.05  0.67  0.7 0.1]);
       % title('channel 3  Data');
        v.hAxesRaw4  = axes('position', [0.05  0.56  0.7 0.1]);
       % title('channel 4  Data');
        v.hAxesRaw5  = axes('position', [0.05  0.45  0.7 0.1]);
       % title('channel 5  Data');
        v.hAxesRaw6  = axes('position', [0.05  0.34  0.7 0.1]);
       % title('channel 6  Data');
        v.hAxesRaw7  = axes('position', [0.05  0.23  0.7 0.1]);
      %  title('channel 7  Data');
        v.hAxesRaw8  = axes('position', [0.05  0.12  0.7 0.1]);
       % title('channel 8  Data');
        
        v.hAxesfft= axes('position', [0.78  0.15  0.2 0.80]);   
        title('Processed Data');
        %  end;
        %启动传输
        v.hStartButton = uicontrol('Style', 'PushButton','String', 'Start Acquisition','units', 'pixels','position', [5 10 100 20],'callback', {@startStopCallback});
        %保存数据
        v.hSaveCheck = uicontrol('Style','checkbox','units', 'pixels','position', [110 10 80 20],'string','Save Data','Value',v.saveDataDefault);
        %频率
        v.hHzText = uicontrol('Style', 'text','units', 'pixels','position', [200 8 20 20],'String', 'Hz:', 'backgroundcol', get(gcf, 'color'));
        %1000 HZ default
        v.hHzPopup = uicontrol('Style', 'popupmenu','units', 'pixels','position', [221 10 80 20],'String', v.hzChoices,'Value' , v.hzDefaultIndex);
        %通道数
        v.hChannelText = uicontrol('Style', 'text','units', 'pixels','position', [300 8 50 20],'String', 'Channels:', 'backgroundcol', get(gcf, 'color'));
        %4通道 default
        v.hChannelPopup = uicontrol('Style', 'popupmenu','units', 'pixels','position', [351 10 70 20],'String', v.chChoices,'Value' , v.chDefaultIndex);
        
        v.hModelText = uicontrol('Style', 'text','units', 'pixels','position', [800 8 100 20],'String', '模式选择:', 'backgroundcol', get(gcf, 'color'));
        v.hModel = uicontrol('Style', 'popupmenu','units', 'pixels','position', [900 8 100 20],'String', ['准备模式'; '测试模式'], 'backgroundcol', get(gcf, 'color'));
        
        if (ispc) %provide list of serial ports
            if (exist('winqueryreg') == 0) %this will always work...
                ser = cellstr(strvcat('Simulate data', 'COM1', 'COM2', 'COM3', 'COM4', 'COM5', 'COM6', 'COM7', 'COM8', 'COM9','COM10','COM11','COM12','COM16'));
                if v.serDefaultIndex == 0
                    v.serDefaultIndex = 1;
                end;
            else %Below... hopefully this also works
                %[ok, ser] = system('powershell [System.IO.Ports.SerialPort]::getportnames()'); %<-good in theory, too slow in practice
                %see http://www.mathworks.com/matlabcentral/newsreader/view_thread/265836
                devices = winqueryreg('name', 'HKEY_LOCAL_MACHINE', 'HARDWARE\DEVICEMAP\SERIALCOMM');
                if length(devices) > 0
                    ser = [];
                    for (d = 1:length(devices))
                        ser = [winqueryreg('HKEY_LOCAL_MACHINE', 'HARDWARE\DEVICEMAP\SERIALCOMM', devices{d}), ser ];
                    end;
                    
                    ser = regexp(ser,'\s','split');
                    ser = ser(~cellfun(@isempty, ser)); %deblank
                    ser = ['Simulate data', ser ];
                else
                    ser = cellstr(strvcat('Simulate data'))
                end %devices not empty
                
                if v.serDefaultIndex == 0
                    v.serDefaultIndex = length(ser);
                end;
            end;
        else
            [ok, ser] = system('ls /dev/cu.*');
            ser = regexp(ser,'\s','split');
            ser = ser(~cellfun(@isempty, ser)); %deblank
            ser = ['Simulate data', ser ];
            if (v.serDefaultIndex == 0) %Arduino/Teensy have names like /dev/cu.usbmodem12341 or /dev/cu.usbmodemfa131, our bluetooth will be named /dev/cu.us922kBT0000
                Index = find(not(cellfun('isempty', strfind(ser, 'us'))));
                if isempty(Index)
                    v.serDefaultIndex = 1;
                else
                    v.serDefaultIndex = Index(1); %choose first port that matches our search string
                end;
            end; %serDefaultIndex = 0
        end; %if ispc else
        if v.serDefaultIndex >  length(ser) %user specified a port that does not exist
            v.serDefaultIndex =  length(ser);
        end;
        v.hSerialPopup = uicontrol('Style', 'popupmenu','units', 'pixels','position', [560 10 200 20],'String', ser,'Value' , v.serDefaultIndex);
        
        % TCP server
        v.tcpserver = tcpip('0.0.0.0', 8080, 'NetworkRole', 'server');
        set(v.tcpserver, 'InputBufferSize', 6553600);
        set(v.hStartButton, 'callback', {@startStopCallback});
    %end %makeGUI()


    %%---------------启动传输---------------------------
    function[OK] =  startDeviceType()
        global v;
        v.gUnsavedData = [];
        v.rawData = []; %raw data from device - data that still needs to be decoded into discrete samples
        %v.newSamples = []; %samples acquired in most recent screen refresh
        %采样率 500HZ default
        v.gOscHz = v.hzChoices(get(v.hHzPopup,'Value'));
        %通道数 4通道 default
        v.gOscChannels = v.chChoices(get(v.hChannelPopup,'Value'));
        %采样间隔
        v.gSecPerSample = 1/v.gOscHz;
        %图中显示的点数(2s*500)
        v.gGraphSamples = round(v.gOscHz * v.gGraphTotalTimeSec);
        %0到2之间的1000个点
        v.xData = linspace(0,v.gGraphSamples/v.gOscHz,v.gGraphSamples);
        v.hAxisRaw = [0 v.xData(end) -1 1];
        v.gSampleNumber = 1;
        if get(v.hSaveCheck,'Value') == 1
            %guoyi 0805
            %v.gSaveDataBaseFilename = saveBrainVisionSub([]);
            fprintf('Saving data as %s\n',v.gSaveDataBaseFilename);
        else
            v.gSaveDataBaseFilename = [];
            fprintf('Warning: data not being saved to disk\n');
        end %if gSaveData
        [OK] =startDeviceType2();
%         if (get(v.hSerialPopup,'Value') == 1)
%             v.deviceType = 0; %simulated dataa
%             [OK] =startDeviceType0();
%         else
%             v.deviceType = 1; %myo
%             list=get(v.hSerialPopup,'String');
%             val=get(v.hSerialPopup,'Value');
%             [OK] =startDeviceType1 (list{val});
%         end;
        %v.hAxisRaw('auto y'); %v.hAxisRaw(3) = 0;  v.hAxisRaw(4) = 0;
        % 4*2000的向量
        v.yData = zeros(v.gOscChannels,v.gGraphSamples);
        
        %v.totalData = zeros(v.gOscChannels-1,1);
    %end %startDeviceType()


    %%--启动模拟数据---
    function [OK] = startDeviceType0() %DeviceType 0 = simulated data
        global v;
        fprintf('Simulated Data: Recording %d channels at %d Hz\n',v.gOscChannels,v.gOscHz);
        v.hAxisRaw(3) = -1.5;  v.hAxisRaw(4) = 1.5; %simulated signal in range -1..+1 plus noise
        OK = true;
    %end %startDeviceType0()
    
     %%--通过tcp发送命令---
    function tcp_send_command(cmd);
        global v;
        pause(0.001);
        fwrite(v.tcpserver,[cmd]);
    %end adc_send_command


    %%--通过串口发送命令---
    function adc_send_command(cmd);
        global v;
        pause(0.001);
        fwrite(v.serialObj,[cmd]);
    %end adc_send_command
    
    function adc_wreg(reg, val)
        %see pages 40,43 of ads1298r datasheet -
        global v;
        pause(0.001);
        fwrite(v.serialObj,[64+reg,0,val]);
    %end adc_wreg


    function [val] =  adc_rreg(reg) %read register
        val = -1;
        global v;
        if (v.serialObj.BytesAvailable > 0)
            fread(v.serialObj, v.serialObj.BytesAvailable); %flush buffer
        end;
        fwrite(v.serialObj,[32+reg,0]);
        for c=1:100%there is often some latency with USB communications
            if (v.serialObj.BytesAvailable < 1)
                pause(0.001);
            end
        end
        if (v.serialObj.BytesAvailable > 0)
            [data,count] = fread(v.serialObj,v.serialObj.BytesAvailable,'uchar');
            val = data(1);
        end
    %end adc_rreg


    function [nMaxChan] = ads_getMaxChannels
        nMaxChan = 8; %ads1294
        fprintf('Channels available: %d  n',nMaxChan);
    %end  ads_getMaxChannels
    
    %%--启动pianote的设备传输----
    function [OK] = startDeviceType2 (DeviceName) %DeviceType 1 =myo on serial port
        %    bandwidth = (gOscHz * ((gOscChannels*3)+1)) /(8 * kBPS); %packets are 4 bytes plus 2 bytes per analog channel  *8 = 8 bits per byte
        global v;
        fopen(v.tcpserver); %等待连接
        %打开连接
        data =[uint8(hex2dec('24')) uint8(hex2dec('24')) uint8(hex2dec('10'))];
        fwrite(v.tcpserver,data);
        OK = true;
        tic;  %计时开始
        
        %同时开始录音   
    %end %startDeviceType2()


    %%--启动myo的设备传输----
    function [OK] = startDeviceType1 (DeviceName) %DeviceType 1 =myo on serial port
        %    bandwidth = (gOscHz * ((gOscChannels*3)+1)) /(8 * kBPS); %packets are 4 bytes plus 2 bytes per analog channel  *8 = 8 bits per byte
        global v;
        RDATAC = uint8(hex2dec('1'));
        OK = false; % assume failure to connect
        if (v.gOscChannels < 1)
            fprintf('Error in %s: set gOscChannels to be at least 1\n',mfilename('fullpath'));
            return;
        end;
        %v.hAxisRaw(3) = -9000000;  v.hAxisRaw(4) = 9000000; %signed 24-bit acquisition data ranges -8388608..8388607
        %开始录音
        if v.audioOn
            record(v.recorder);
        end;
        v.serialObj=serDeviceIndexSub (DeviceName);
        maxChan = ads_getMaxChannels;
        if (maxChan < 1)
            fprintf('Error in %s: unable to detect an ads1298-compatible system. Please check connections.\n',mfilename('fullpath'));
            return;
        end;
        if (maxChan < v.gOscChannels)
            fprintf('Error ads1298-compatible system can only support up to %d channels.\n',maxChan);
            return;
        end;
        OK = true;
        if isempty(DeviceName)
            fprintf('ads Data: Recording %d channels at %d Hz\n',v.gOscChannels-1,v.gOscHz);
        else
            fprintf('ads Data: Recording %d channels at %d Hz attached to port "%s"\n',v.gOscChannels-1,v.gOscHz,DeviceName);
        end;
        if (v.testSignal)
            adc_send_command(TEST_SIGNAL);
        end;
        v.gOscChannels = v.gOscChannels+1; %include extra channel for digital data    5  guoyi  现在没用
        
        adc_send_command(RDATAC); %restart streaming
        tic;  %计时开始
        
        %同时开始录音
        
        
    %end %startDeviceType1()


    %%-停止设备传输-----------------------
    function stopDeviceType()
        global v;
        stopDeviceType2();
%         switch v.deviceType
%             case 1
%                 stopDeviceType1();
%             otherwise
%                 stopDeviceType0();
%         end %switch deviceType
        flushSaveDataSub(); %save any residual data to disk
        fcloseSerialSub();
        totalanalysis();
    %end %stopDeviceType()


    function stopDeviceType0(v)
        disp('Simulated acquisition halted');
    %end %stopDeviceType0()
    
    function stopDeviceType2()
        global v;
        disp('pianote acquisition halted');
        data =[uint8(hex2dec('24')) uint8(hex2dec('24')) uint8(hex2dec('00'))];
        fwrite(v.tcpserver,data);
        fprintf('肌肉电记录总时间 %ss  \n',num2str(toc));
        fprintf('记录的数据总长度: mpu:%s s  emg: %s s\n',num2str(length(v.mpuData)/100),num2str(length(v.totalData)/1000));


        pause(1);
%         data =uint8(hex2dec('23'));
%         fwrite(v.tcpserver,data);
        pause(1);
      
    %%---------------------------------------------------
    function stopDeviceType1()
        global v;
        disp('myo acquisition halted');
        if ~isnumeric(v.serialObj) &&  isvalid(v.serialObj)
            SDATAC = uint8(hex2dec('2'));
            adc_send_command(SDATAC); %stop any sampling - we will configure setup
            fprintf('肌肉电记录总时间 %ss  \n',num2str(toc));
            %同步停止录音,将文件存储
            if v.audioOn
                stop(v.recorder);
                load chirp.mat;
                v.myRecording = getaudiodata(v.recorder);
                fprintf('声音记录总时间 %ss  \n',num2str(length(v.myRecording)/44100));
                v.datestr = datestr(now,'yymmdd_HHMMSS');     
                filename = [v.saveaudiofile v.datestr '.wav'];
                audiowrite(filename,v.myRecording,v.audiofs);
            end;    
        end;
    %end %stopDeviceType1()


    %%---------------------------------------------------
    function flushSaveDataSub()
        global v;
        if isempty(v.gSaveDataBaseFilename)
            return;
        end
        v.gSaveNumber = 1;
       % saveBrainVisionSub(v.gUnsavedData,v.gOscHz,false, true, v.gSaveDataBaseFilename);
        v.gUnsavedData = [];
    %end; %flushSaveDataSub()






    %%----------启动传输的处理函数-----------------------------
    function startStopCallback(hObject, eventdata)
        global v;
        if v.acquiringData
            if strcmp(v.timerObj.running, 'on')
                stop(v.timerObj);
            end
            stopDeviceType;
            v.acquiringData = false;
            set(hObject, 'string', 'Start Acquisition');
            controlStatus = 'on';
        else
            OK = startDeviceType();
            if OK == false
                return;
            end;
            v.acquiringData = true;
            set(hObject, 'string', 'Stop Acquisition');
            controlStatus = 'off';
            if strcmp(v.timerObj.running, 'off')
                start(v.timerObj);
            end
        end
        set(v.hSaveCheck,'Enable',controlStatus);
        set(v.hSerialPopup,'Enable',controlStatus);
        set(v.hChannelPopup,'Enable',controlStatus);
        set(v.hHzPopup,'Enable',controlStatus);
    %end %startStopCallback()


    %%---------------------------------------------------
    function figureCloseCallback(hObject, eventdata)
        global v;
        cleanupObjects();
    %end %figureCloseCallback()


    %%---------------------------------------------------
    function cleanupObjects()
        global v;
        if isvalid(v.timerObj)
            stop(v.timerObj);
            delete(v.timerObj);
        end
        fcloseSerialSub();
        if ishandle(v.hFigure),
            delete(v.hFigure);
        end
    %end %cleanupObjects()


    %%---------------------------------------------------
    function fcloseSerialSub()
        global v;
        if ~isnumeric(v.serialObj) &&  isvalid(v.serialObj)
            if (v.serialObj.BytesAvailable > 0)
                fread(v.serialObj, v.serialObj.BytesAvailable); %flush buffer
            end;
          %  SDATAC = uint8(hex2dec('2'));
          %  adc_send_command(SDATAC); %stop any sampling - we will configure setup
            
            fclose(v.serialObj);
            delete(v.serialObj);
            disp('closed serial port');
        end
        fclose(v.tcpserver);
    %end %fcloseSerialSub()


    %%---------------------------------------------------
    function [theSamples, rawResidual] =serDecodeSub(rawData)
        global v;
        theSamples =[];
        mpuSample = [];
        qSample = []; %四元数
        euSample = [];
        accSample = [];
        len = length(rawData);
        samples = 0;
        % 0x24|type|data
        %|type标示这个包中包含的数据,有ads1294和mpu6000两种。暂时只分析了ads1294一种 
        pos = 1;
        OKpos = 0;
        while ((len-pos+1) >= 64)  %仅分析EMG数据 其他类型的比这长度小的数据包也放在下次分析。
            %new protocol guoyi .10.10 
           % fprintf('Len %d rawData %d: %d \n',len , pos,rawData(pos) );
            if (rawData(pos) == 36 ) %header byte always starts 0x24
              % fprintf('rawData(pos)  %d \n',rawData(pos) );
               pos = pos + 1;
               if (rawData(pos) == 36 ) %header
                    pos = pos + 1;
                    if(rawData(pos) == 1) % DATA_BIO_GEST
                        pos =pos +1; %长度,暂时不检查,guoyi. 2014-11-10
                       % fprintf('gesture %d\n',rawData(pos+1));
                       % if(rawData(pos+1) ~=6)
                       %     set(v.hGestureText,'String',rawData(pos+1));
                       % end;
                        pos = pos +1 +1;
                        OKpos = pos -1;
                      %  continue;
                    end;
                    if(rawData(pos) == 2) % DATA_MPU_RAW ,ax,ay,az,gx,gy,gz
                        %fprintf('DATA_MPU_IMU\n');
                        pos =pos +1; %长度,暂时不检查,guoyi. 2014-11-10
                        mpuSample(1,1) = typecast(uint8([rawData(pos+2) rawData(pos+1)]),'int16');
                        mpuSample(2,1) = typecast(uint8([rawData(pos+4) rawData(pos+3)]),'int16');
                        mpuSample(3,1) = typecast(uint8([rawData(pos+6) rawData(pos+5)]),'int16');
                        mpuSample(4,1) = typecast(uint8([rawData(pos+8) rawData(pos+7)]),'int16');
                        mpuSample(5,1) = typecast(uint8([rawData(pos+10) rawData(pos+9)]),'int16');
                        mpuSample(6,1) = typecast(uint8([rawData(pos+12) rawData(pos+11)]),'int16');
                        v.mpuData =[v.mpuData mpuSample];
                        pos  = pos + 2*6 +1; % ignore 9轴原始数据. 
                        OKpos = pos -1;
                     %   continue;   %next
                    end;
                    if(rawData(pos) == 3) % DATA_MPU_EU
                       % fprintf('DATA_MPU_EU\n');
                        pos =pos +1; %长度,暂时不检查,guoyi. 2014-11-10
                        euSample(1,1) = typecast(uint8([rawData(pos+5) rawData(pos+4) rawData(pos+3) rawData(pos+2)]),'single');
                        euSample(2,1) = typecast(uint8([rawData(pos+9) rawData(pos+8) rawData(pos+7) rawData(pos+6)]),'single');
                        euSample(3,1)  = typecast(uint8([rawData(pos+13) rawData(pos+12) rawData(pos+11) rawData(pos+10)]),'single');
                        v.euData =[v.euData euSample];
                        pos  = pos + 4*3 +2; % ignore 惯导欧拉角. 
                        OKpos = pos -1;
                    %    continue;   %next
                    end;
                    if(rawData(pos) == 4) % DATA_MPU_ACC
                       % fprintf('DATA_MPU_ACC\n');
                        pos =pos +1; %长度,暂时不检查,guoyi. 2014-11-10
                        accSample(1,1) = typecast(uint8([rawData(pos+3) rawData(pos+2)]),'int16');
                        accSample(2,1) = typecast(uint8([rawData(pos+5) rawData(pos+4)]),'int16');
                        accSample(3,1) = typecast(uint8([rawData(pos+7) rawData(pos+6)]),'int16');
                         
                        v.accData =[v.accData accSample];
                        pos  = pos + 3*2 +2; % ignore 3轴向加速度.   
                        OKpos = pos -1;
                   %     continue;   %next
                    end;
                    if(rawData(pos) == 5) % 四元数
                        %fprintf('DATA_MPU_EU\n');
                        pos =pos +1; %长度,暂时不检查,guoyi. 2014-11-10
                        qSample(1,1) = typecast(uint8([rawData(pos+4) rawData(pos+3) rawData(pos+2) rawData(pos+1)]),'int32');
                        qSample(2,1) = typecast(uint8([rawData(pos+8) rawData(pos+7) rawData(pos+6) rawData(pos+5)]),'int32');
                        qSample(3,1)  = typecast(uint8([rawData(pos+12) rawData(pos+11) rawData(pos+10) rawData(pos+9)]),'int32');
                        qSample(4,1)  = typecast(uint8([rawData(pos+16) rawData(pos+15) rawData(pos+14) rawData(pos+13)]),'int32');
                        v.qData =[v.qData qSample];
                        pos  = pos + 4*4 +1; % ignore 惯导欧拉角. 
                        OKpos = pos -1;
                    %    continue;   %next
                    end;
                    if(rawData(pos) == 8) % DATA_EMG_RAW
                        %fprintf('DATA_EMG_RAW\n');
                        pos =pos +1; %长度,暂时不检查,guoyi. 2014-11-10
                        samples = samples + 1; %采样值++
                        %fprintf('channel: %d samples %d: ',v.gOscChannels, samples);
                        for i=1:6
                             theSamples(i,samples) = typecast(uint8([rawData(pos+(i*2)) rawData(pos-1+(i*2))]), 'uint16'); %matlab does not natively support signed 24-bit, so create signed 32-bit
                           % fprintf('%x  %x %x \n',rawData(pos-2+(i*3)), rawData(pos-1+(i*3)), rawData(pos+(i*3)));
                            %theSamples(i,samples) = typecast(uint8([rawData(pos+(i*3)) rawData(pos-1+(i*3)) rawData(pos-2+(i*3))  0]), 'uint32'); %matlab does not natively support signed 24-bit, so create signed 32-bit


                            %theSamples(i,samples) = bitshift(theSamples(i,samples),-8,'int32'); %convert 32-bit to 24-bit
                            %fprintf('%d  \n',theSamples(i,samples));
                        end;
                        %fprintf('\n');
                       % theSamples(v.gOscChannels,samples) = bitand(rawData(pos),15);
                        pos = pos + 2*6 + 1; % EMG 8channel data
                        OKpos = pos -1;%last valid byte
                    end
                    if(rawData(pos) == 16) % DATA_EMG_FILT , 进过嵌入式系统内过滤后的EMG数据
                        pos =pos +1; %长度,暂时不检查,guoyi. 2014-11-10
                        samples = samples + 1; %采样值++
                        %  fprintf('samples %d: ',samples);
                        for i=1:(v.gOscChannels-1)
                            theSamples(i,samples) = typecast(uint8([rawData(pos+(i*4)) rawData(pos-1+(i*4)) rawData(pos-2+(i*4)) rawData(pos-3+(i*4))]),'int32');
                            %theSamples(i,samples) = typecast(uint8([0 rawData(pos+(i*3)) rawData(pos-1+(i*3)) rawData(pos-2+(i*3)) ]), 'int32'); %matlab does not natively support signed 24-bit, so create signed 32-bit
                            %theSamples(i,samples) = bitshift(theSamples(i,samples),-8,'int32'); %convert 32-bit to 24-bit
                            %fprintf('%x  ',theSamples(i,samples));
                        end;
                         %fprintf('\n');
                        theSamples(v.gOscChannels,samples) = bitand(rawData(pos),15);
                        pos = pos + 4*v.gOscChannels + 1; % EMG 8channel data
                        OKpos = pos -1;%last valid byte
                    end
                else
                    pos = pos + 1;
                end; %if command else data
            else
                 pos = pos + 1;
            end;
        end;
        if (OKpos < len)
            %We need to retain left over bytes. Consider if the serial port transfers blocks of 512 bytes, but the expected packet
            %sise is 6 bytes, there will always be left over bytes
            %fprintf(' residual bytes: decoded %d of %d bytes (packet size %d)\n', OKpos,len,packetBytes );
            rawResidual = rawData((OKpos+1):end);
        else
            rawResidual = [];
        end;
        %datestr(now)
        v.totalsamples = v.totalsamples + samples;
        fprintf('Decoding:  add %d simples\n',samples);
    %end; % serDecodeSub()


    %%---------------------------------------------------
    function ser=serDeviceIndexSub(DeviceName)
        if ~nargin || isempty(DeviceName)
            if (ispc)
                DeviceName = 'COM9';
            else
                %DeviceName = '/dev/tty.usbmodem12341';
                DeviceName = '/dev/cu.usbmodem12341';
            end;
        end; %DeviceName not specified
        if (ispc)
            fprintf('Assuming device is named "%s", use the Device Manager to show active ports.\n', DeviceName);
        else
            fprintf('Assuming device is named "%s", available port names are\n', DeviceName);
            system('ls /dev/cu.*');
        end;
        %ser = serial(DeviceName,'InputBufferSize',1638400,'BaudRate',921600);%most reliable to set correct BPS - helps if device is already running at high speeds
        ser = serial(DeviceName,'InputBufferSize',1638400,'BaudRate',38400);
        fopen(ser);
    %end % SerDeviceIndexSub()
    
    function  totalanalysis( )
        %ANALYSIS 此处显示有关此函数的摘要
        %   此处显示详细说明
        global v;
        
        if  length(v.totalData) < 3000
            return;
        end;
        data = double(v.totalData);
        data = data';
        if v.testmodel 
            v.datestr = '';
        end
        filename = [v.saverawfile v.datestr 'raw.txt'];
        save(filename,'data', '-ASCII');
        if isempty(v.euData) == 0
            fprintf('end...\n');
            eudata = v.euData';
            filename = [v.savefiltfile v.datestr 'eu.txt'];
            save(filename,'eudata', '-ASCII');
        end
        if isempty(v.accData) == 0
            accdata = v.accData';
            filename = [v.savefiltfile v.datestr 'acc.txt'];
            save(filename,'accdata', '-ASCII');
        end
        if isempty(v.qData) == 0
            qdata = v.qData';
            filename = [v.savefiltfile v.datestr 'quat.txt'];
            save(filename,'qdata', '-ASCII');
        end
        if isempty(v.mpuData) == 0
            mpudata = v.mpuData';
            filename = [v.savefiltfile v.datestr 'mpu.txt'];
            save(filename,'mpudata', '-ASCII');
        end
        fprintf('记录的数据总长度: mpu:%s s  emg: %s s\n',num2str(length(v.mpuData)/100),num2str(length(data)/1000));
        
        data = data';
        data = data/4096*3.3;
        %data = data*3.3/4096;
        data = filtfilt(v.b,1,data');
        data = filtfilt(v.b2,1,data);
        data = filtfilt(v.b1,1,data);
        data = data(:,v.channel_useful)';
        rawdata = data';
        filename = [v.savefiltfile v.datestr 'filt.txt'];
        save(filename,'rawdata', '-ASCII');
        
        %分析
        if v.testmodel == false
            fprintf('启动弹奏分析....\n');
            flag =  get(v.hModel,'Value');
            main(qdata,mpudata,v.myRecording,rawdata,flag);
        end
                
        v.totalData = [];
        v.totalsamples  = 0;
        v.accData = [];
        v.euData = [];
        v.qData = [];
        v.mpuData = [];
        if v.testmodel 
            testdata;
        end


        
%         figure(2);
%         set(gcf,'outerposition',get(0,'screensize'));
%         subplot(8,1,1);
%         plot(data(1,:));
%         title('filter data using matlab');
%         subplot(8,1,2);
%         plot(data(2,:));
%         subplot(8,1,3);
%         plot(data(3,:));
%         subplot(8,1,4);
%         plot(data(4,:));
%         subplot(8,1,5);
%         plot(data(5,:));
%         subplot(8,1,6);
%         plot(data(6,:));
%         subplot(8,1,7);
%         plot(data(7,:));
%         subplot(8,1,8);
%         plot(data(8,:));
%         figure(3);
%         [freq,fftdata] = powerSpectrum(length(data(1,:)), data(v.showPowerChannel,:));
%         plot(freq, fftdata);
%         xlabel('Frequency (Hz)'); ylabel('Amplitude');
        %save('data.mat','data');


        
        %analysis(data);
    %%end
    

展开阅读全文

没有更多推荐了,返回首页