与傅里叶的故事

什么都不说，先放一个肌电手环的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.gSe