Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Python into matlab

Subject: Python into matlab

From: Armindo

Date: 11 May, 2014 03:45:13

Message: 1 of 1

Please can someone help me to convert this python code into matlab code.
Any help would be apreciated: The code come from this page: http://nbviewer.ipython.org/github/duartexyz/BMC/blob/master/ResidualAnalysis.ipynb

"""Automatic search of filter cutoff frequency based on residual analysis."""

from __future__ import division, print_function
import numpy as np
from scipy.signal import butter, filtfilt

def residual_analysis(y, freq=1, fclim=[], show=False):

from scipy.interpolate import UnivariateSpline

    # signal filtering
    freqs = np.linspace((freq / 2) / 101, freq / 2, 101)
    res = []
    for fc in freqs:
        b, a = butter(2, fc / (freq / 2))
        yf = filtfilt(b, a, y)
        # residual between filtered and unfiltered signals
        res = np.hstack((res, np.sqrt(np.mean((yf - y) ** 2))))

    # find the optimal cutoff frequency by fitting an exponential curve
    # y = A*exp(B*x)+C to the residual data and consider that the tail part
    # of the exponential (which should be the noisy part of the residuals)
    # decay starts after 3 lifetimes (exp(-3), 95% drop)
    if not len(fclim) or np.any(fclim < 0) or np.any(fclim > freq / 2):
        fc1 = 0
        fc2 = np.nonzero(freqs >= 0.9 * (freq / 2))[0][0]
        # log of exponential turns the problem to first order polynomial fit
        # make the data always greater than zero before taking the logarithm
        reslog = np.log(np.abs(res[fc1:fc2 + 1] - res[fc2]) +
                        10 * np.finfo(np.float).eps)
        Blog, Alog = np.polyfit(freqs[fc1:fc2 + 1], reslog, 1)
        fcini = np.nonzero(freqs >= -3 / Blog) # 3 lifetimes
        fclim = [fcini[0][0], fc2] if np.size(fcini) else []
    else:
        fclim = [np.nonzero(freqs >= fclim[0])[0][0],
                 np.nonzero(freqs >= fclim[1])[0][0]]

    # find fc_opt with linear fit y=A+Bx of the noisy part of the residuals
    if len(fclim) and fclim[0] < fclim[1]:
        B, A = np.polyfit(freqs[fclim[0]:fclim[1]], res[fclim[0]:fclim[1]], 1)
        # optimal cutoff frequency is the frequency where y[fc_opt] = A
        roots = UnivariateSpline(freqs, res - A, s=0).roots()
        fc_opt = roots[0] if len(roots) else None
    else:
        fc_opt = None

    if show:
        _plot(y, freq, freqs, res, fclim, fc_opt, B, A)

    return fc_opt


def _plot(y, freq, freqs, res, fclim, fc_opt, B, A):
    """Plot results of the residual_analysis function, see its help."""
    try:
        import matplotlib.pyplot as plt
    except ImportError:
        print('matplotlib is not available.')
    else:
        #fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(6, 8))
        plt.figure(num=None, figsize=(10, 5))
        ax1 = plt.subplot(121)
        plt.rc('axes', labelsize=12, titlesize=12)
        plt.rc('xtick', labelsize=12)
        plt.rc('ytick', labelsize=12)
        ax1.plot(freqs, res, 'b.', markersize=9)
        time = np.linspace(0, len(y) / freq, len(y))
        ax2 = plt.subplot(222)
        ax2.plot(time, y, 'g', linewidth=1, label='Unfiltered')
        ydd = np.diff(y, n=2) * freq ** 2
        ax3 = plt.subplot(224)
        ax3.plot(time[:-2], ydd, 'g', linewidth=1, label='Unfiltered')
        if fc_opt:
            ylin = np.poly1d([B, A])(freqs)
            ax1.plot(freqs, ylin, 'r--', linewidth=2)
            ax1.plot(freqs[fclim[0]], res[fclim[0]], 'r>',
                     freqs[fclim[1]], res[fclim[1]], 'r<', ms=10)
            ax1.set_ylim(ymin=0, ymax=4 * A)
            ax1.plot([0, freqs[-1]], [A, A], 'r-', linewidth=2)
            ax1.plot([fc_opt, fc_opt], [0, A], 'r-', linewidth=2)
            ax1.plot(fc_opt, 0, 'ro', markersize=7, clip_on=False,
                     zorder=9, label='$Fc_{opt}$ = %.1f Hz' % fc_opt)
            ax1.legend(fontsize=12, loc='best', numpoints=1, framealpha=.5)
            b, a = butter(2, fc_opt / (freq / 2))
            yf = filtfilt(b, a, y)
            ax2.plot(time, yf, color=[1, 0, 0, .5],
                     linewidth=2, label='Opt. filtered')
            ax2.legend(fontsize=12, loc='best', framealpha=.5)
            ax2.set_title('Signals (RMSE = %.3g)' % A)
            yfdd = np.diff(yf, n=2) * freq ** 2
            ax3.plot(time[:-2], yfdd, color=[1, 0, 0, .5],
                     linewidth=2, label='Opt. filtered')
            ax3.legend(fontsize=12, loc='best', framealpha=.5)
            resdd = np.sqrt(np.mean((yfdd - ydd) ** 2))
            ax3.set_title('Second derivatives (RMSE = %.3g)' % resdd)
        else:
            ax1.text(.5, .5, 'Unable to find optimal cutoff frequency',
                     horizontalalignment='center', color='r', zorder=9,
                     transform=ax1.transAxes, fontsize=12)
            ax2.set_title('Signal')
            ax3.set_title('Second derivative')

        ax1.set_xlabel('Cutoff frequency [Hz]')
        ax1.set_ylabel('Residual RMSE')
        ax1.set_title('Residual analysis')
        ax1.grid()
        #ax2.set_xlabel('Time [s]')
        ax2.set_xlim(0, time[-1])
        ax2.grid()
        ax3.set_xlabel('Time [s]')
        ax3.set_xlim(0, time[-1])
        ax3.grid()
        plt.tight_layout()
        plt.show()

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us