% settings
SerialPort='/dev/tty.usbmodem641'; % serial port
N = 500;
Fs = 100;

m=zeros(1,N);

% initialize
figure(1)
hLine = plot(m);
ylim([0 1024]);
set(hLine,'YData',m);

figure(2)
L=length(m);
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(m,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);

% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1))) 
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')

KeepRunning = 1;
while KeepRunning
    s = serial(SerialPort);
    set(s,'BaudRate',57600);
    fopen(s);

    for i = 1:N
        datum = fscanf(s, '%s');
        fprintf('%s\n', datum);

        if (length(datum) > 0)
            m(i) = str2num(datum);
        else
            m(i) = 0;
        end
    end

    % Clean up the serial port
    fclose(s);
    delete(s);
    clear s;
    
    f0 = 87.5;                %#notch frequency
    fn = Fs/2;              %#Nyquist frequency
    freqRatio = f0/fn;      %#ratio of notch freq. to Nyquist freq.

    notchWidth = 1;       %#width of the notch

    %Compute zeros
    zeros = [exp( sqrt(-1)*pi*freqRatio ), exp( -sqrt(-1)*pi*freqRatio )];

    %#Compute poles
    poles = (1-notchWidth) * zeros;

    figure(3);
    zplane(zeros.', poles.');

    b = poly( zeros ); %# Get moving average filter coefficients
    a = poly( poles ); %# Get autoregressive filter coefficients
 
    figure(4);
    freqz(b,a,32000,Fs)

    %#filter signal x
    y = filter(b,a,m);

    %y = m;
    
    % Remove DC offset
    mu = mean(y);
    y = y - mu;

    figure(1)
    hLine = plot(y(1:end));
    ylim([-512 512]);
    set(hLine,'YData',y);

    figure(2)
    L=length(y);
    NFFT = 2^nextpow2(L); % Next power of 2 from length of y
    Y = fft(y,NFFT)/L;
    f = Fs/2*linspace(0,1,NFFT/2+1);

    % Plot single-sided amplitude spectrum.
    plot(f,2*abs(Y(1:NFFT/2+1))) 
    title('Single-Sided Amplitude Spectrum of y(t)')
    xlabel('Frequency (Hz)')
    ylabel('|Y(f)|')
end