Testing fsolve and anonymous function

I have got a need to solve equations repeatedly whose coefficients are changing with each iteration. MatLab function 'fsolve' is tested using function handle and anonymous function as follows.

s=[3 7 11];
d=[1 3 5];
%ans= (2,1), (5,2), (8,3)
for i=1:3
x0 = [-5; -5]; % Make a starting guess at the solution
myfun =@(x) [x(1) + x(2) - s(i);
x(1) - x(2) - d(i)];
[x,fval] = fsolve(myfun,x0)  % Call 
end

Hand, Foot and Mouth Disease (HFMD)

My baby got blister-like rash on her hands and feet last week. Doctor said it may be caused by HFMD or flu. I found out later that it was the flu in my daughter case. I have never heard of HFMD before and that is why I re-share a brochure from KK hospital.

Q. What is hand, foot and mouth disease (HFMD)?
A. This is an infectious disease caused by a family of viruses called Enteroviruses, the commonest being the Coxsackie virus and Enterovirus. It can occur in people from various age groups, especially in pre-schoolers. It is a very common disease in Singapore and has been in existence for many years. It is not a rare or new disease.

Q. How do you know if your child has HFMD?
A. Children with HFMD will have blister-like rash on hands, feet and buttocks, mouth ulcers and fever. In addition, the child may have a sore throat, runny nose, vomiting and diarrhoea, and may feel tired. You may bring your child to the polyclinic or see your family doctor. There is no need for you to rush your child to the Children's Emergency just to confirm the diagnosis of HFMD.

Q. How can your child get HFMD?
A. HFMD can be easily spread through direct contact with nose discharge, saliva, faeces and fluid from the blisters.

Q. Is this disease serious?
A. The disease is usually mild and most children will recover in about a week's time. Only very rarely do certain rare strains of the virus cause complications such as inflammation of the brain and heart.

Q. Is HFMD treatable?
A. There is no specific treatment for HFMD. The symptoms are usually mild and children usually recover well as their own immune system fights off the virus. Your doctor will give medication to control the fever. You should encourage your child to take as much oral fluids as possible. Your child may not have a good appetite because swallowing may be painful. However, ensure that your child has adequate fluid to prevent dehydration. Offer your child small amount of fluid such as diluted fruit juice, rice or barley water every half hourly and about 10 to 30 ml each time throughout the day. Antibiotics are ineffective because this is a viral, not a bacterial infection.

Q. Does a prior infection with enterovirus make a person immune?
A. Specific immunity can occur, but a second episode is possible from a different strain of virus belonging to the enterovirus family.

Q. What can be done to prevent the spread of this disease?
A. Infected children should not be allowed to go to school, childcare centres and other crowded places until he is fully recovered. Practise good general hygiene. Wash your hands immediately after contact with the infected child or handling diaper changes, and before handling food. Prevent other children from contact with toys, books, eating utensils, towels, clothes and other personal items used by the infected child.

Q. When should a child with HFMD be brought to the Children's Emergency?
A. Most children with HFMD are relatively well and active despite their illness. Your family doctor or the polyclinic will be able to manage the majority of the cases. However, you should bring your child to the Children's Emergency if he develops any of the following problems:
*When the oral intake of fluids is poor, or when the child is unable to swallow, or vomits persistently.
*When the tongue is dry, or when the child has decreased urine output (dehydration).
*If the child appears lethargic, drowsy or irritable, is crying persistently, or is disoriented.
*When seizures occur.
*If there is difficulty in breathing.
*If the child looks ashen, pale or blue.
*If the child complains of acute headache or giddiness, or if there is neck stiffness.

Geometric Template Matching in LabVIEW

In NI's IMAQ Vision Concepts Manual, geometric template matching is described as follows. Geometric matching locates regions in a grayscale image that match a model, or template, of a reference pattern. Geometric matching is specialized to locate templates that are characterized by distinct geometric or shape information. When using geometric matching, a template is created that represents the object to be searched. Machine vision application then searches for instances of the template in each inspection image and calculates a score for each match. The score relates how closely the template resembles the located matches. Geometric matching finds template matches regardless of lighting variation, blur, noise, occlusion, and geometric transformations such as shifting, rotation, or scaling of the template.
The VIs such as IMAQ Find CoordSys (Pattern) 2 are used to locate the model. The template for the model is created as discussed in the following steps.

Open Template Editor in Windows by clicking Start -> All Programs -> National Instruments -> Vision -> Template Editor. Click File menu->New Template.... Select Geometric Matching Template (Edge Based) and browse an image to extract the template from. In the Select Template Region tab, define a region. For example, start at (200,300) and drag the mouse cursor to (232,332) and release it. Then, you can move the selection to the desired location. In the Define curves tab, specify curve parameters, e.g., Extraction mode to normal, Edge Threshold to 32, Edge Filter size to Fine, Minimum Length to 5, Row Search Step Size to 1 and Column Search Step size to 1. The setting in Customize Scoring and Specify Match Options tabs are set as default. Save the template by clicking File -> Save Template....

As an example, I have created a few VIs at the following link.

Geometric Template Matching on GitHub

Initialization for the VI inputs is shown below.

Troubleshooting Daikin Aircon

We have been using Daikin(3MKD58DVM) inverter multi-split (R-22) aircon for a few years. Last month, one of the indoor units (FTKD25DVM) was faulty. There was power and we could turn it on but the operation stops suddenly after a few seconds. Then one the LED indicators kept blinking. After searching on the Internet, we found a troubleshooting procedure using our ARC433B47 remote controller as follows.

Check Method 1

1. When the timer cancel button is held down for 5 seconds, 00 indication appears on the temperature display section.

2. Press the timer cancel button repeatedly until a long beep sounds. A short beep and two consecutive beeps indicate non-corresponding codes. The code indication changes in the sequence shown below.

3. When you hear the long beep, see the code on the indicator. It is the error code and description of fault is shown in the table below.

4. To return to the normal mode, hold the timer cancel button down for 5 seconds. When the remote controller is left untouched for 60 seconds, it also returns to the normal mode.

Check Method 2

1. Press the center of the TEMP button and the MODE button at the same time. The left-side digit will blink.

2. Press the TEMP up or down button and change the figure until you hear the beep.

3. Diagnose by the sound.

beep : The left-side digit does not correspond with the error code.
beep beep : The left digit corresponds with the error code but the right-side digit does not.
long beep : Both digits correspond with the error code.


4. Press the MODE button. The right-side digit will blink.

5. Press the TEMP up or down button and change the figure until you hear the beep.

6. Diagnose by the sound.

beep : The left-side digit does not correspond with the error code.
beep beep : The left digit corresponds with the error code but the right-side digit does not.
long beep : Both digits correspond with the error code.


7. When you hear the long beep, see the code on the indicator. It is the error code and description of fault is shown in the table below.

8. Press the MODE button to exit from the diagnosis mode.

9. Press the ON/OFF button twice to return to the normal mode. When the remote controller is left untouched for 60 seconds, it also returns to the normal mode.

Error Codes and Description of Fault

System

00 Normal
U0 Refrigerant shortage
U2 Low-voltage detection or over-voltage detection
U4 Signal transmission error (between indoor and outdoor units)
U7 Signal transmission error (on outdoor unit PCB)
UA Unspecified voltage (between indoor and outdoor units)
UH Anti-icing function in other rooms

Indoor Unit

A1 Indoor unit PCB abnormality
A5 Freeze-up protection control or heating peak-cut control
A6 Fan motor or related abnormality
C4 Indoor heat exchanger thermistor or related abnormality
C9 Room temperature thermistor or related abnormality

Outdoor Unit

A5 Anti-icing function
E1 Outdoor unit PCB abnormality
E5 OL activation (compressor overload)
E6 Compressor lock
E7 DC fan lock
E8 Input overcurrent detection
EA Four-way valve abnormality
F3 Discharge pipe temperature control
F6 High pressure control in cooling
H0 Compressor sensor system abnormality
H6 Position sensor abnormality
H8 CT or related abnormality
H9 Outdoor temperature thermistor or related abnormality
J3 Discharge pipe thermistor or related abnormality
J6 Outdoor heat exchanger thermistor or related abnormality
J8 Liquid pipe thermistor or related abnormality
J9 Gas pipe thermistor or related abnormality
L3 Electrical box temperature rise
L4 Radiation fin temperature rise
L5 Output overcurrent detection
P4 Radiation fin thermistor or related abnormality

In our case, remote controller display was U4. It means Signal Transmission Error (between Indoor and Outdoor Units). The data sent from the outdoor unit cannot be received normally, or the content of the data is abnormal. Supposed causes are

1. Defective outdoor unit PCB
2. Defective indoor unit PCB
3. Wiring error
4. Disturbed power supply waveform
5. Breaking of the connection wires between the indoor and outdoor units (wire No. 3)


When the service guys came and checked our aircon, it was found that wire 3 was broken due to erosion at the connection point on outdoor unit. They suggested to replace the whole wire set in our system. I thought it was not necessary. After they peeled away one inch of outer layer and connected back, the system was OK. Interestingly, only third wire for all indoor units were eroded. It looks there must be DC voltage in that wire and it may be acting as electrode. So, I think the same thing will happen again.

2013 Jun 10 Update :

The system has been OK so far. It might not be electrode problem according to a comment on this blog post which said Daikin does not use DC for the signal. It might be because of the bad quality wire that the contractor used. I can now guess why they insisted to change the whole wire set :P Thanks to those who wrote comments in this post.

2014 Feb 10 Update :

The same problem happened again in another indoor unit. This time, it was caused by faulty outdoor unit PCB. The said the cost to change the PCB is around 800 SGD.

Reference:
Malfunction code chart
Service Manual

Curve Fitting in Matlab

One of my friends was doing analysis of some sampled data. He wanted to fit them into a curve to derive the equation. He said Excel could fit up to polynomial of order 6 and he wanted a higher order one. So, he asked my help to write a Matlab program. He also wanted the plotted curve besides the calculated coefficients. The program and sample data are shown below. That is good for me because I also need to fit data occasionally and I can reuse it later :)
See d.txt and cf.m.

%--------------------------------------
%Curve fitting
%2011-Jun-03
%Programmer: Yan Naing Aye
%--------------------------------------
clc;
close all;
clear all;
%--------------------------------------
%get data from file that can be edited
load('d.txt');
X=d(:,1);
Y=d(:,2);
%--------------------------------------
%define order
N=2;
%--------------------------------------
%curve fitting
p=polyfit(X,Y,N);
%generate curve fitted data
[nr nc]=size(X);
YA=zeros(nr,1);
for i=0:N
    YA=YA+p(i+1)*X.^(N-i);
end
%--------------------------------------
%Plot Y vs X
figure;
plot(X,Y,' rd','LineWidth',2,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',6)
hold on;
plot(X,YA,'b','LineWidth',2,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor','g',...
                'MarkerSize',6)
hold off;            
grid on;
Title('Y vs X')
xlabel('X');
ylabel('Y');
legend('Sampled data','Fitted curve',...
       'Location','NW')
%--------------------------------------

Monte Carlo integration

As I am reading about Particle Filter, just for fun, I wrote a Matlab program that performs Monte Carlo integration using uniform distribution. The equation and the program are as follows.

clc;
close all;
clear all;
%--------------------------------------
%Monte Carlo Integration
%--------------------------------------
%--------------------------------------
%limits
LowerLimit=1;
UpperLimit=3;
%--------------------------------------
%generate uniformly distributed random 
%numbers within the range
n=10000;
x=LowerLimit+(UpperLimit-LowerLimit)*rand(n,1);
%--------------------------------------
%Probability distribution function for 
%uniformly distributed random numbers
q=1/(UpperLimit-LowerLimit);
%--------------------------------------
%function to integrate
f=3.*x.*x;
i=mean(f./q)
%--------------------------------------

Degrees of Freedom of the Human Arm

Have you ever think how many degrees of freedom the human arm has (excluding palm and fingers)? I would like to share a short excerpt from a book by Saeed Benjamin Niku, Introduction to Robotics -Analysis, Control, Applications. The shoulder, the elbow, and the wrist are three joint clusters in the human arm. The shoulder has 3 degrees of freedom. The upper arm can swing up and down in the coronal plane. It can also swing back and forth in the transverse plane. And rotation with respective to the axis along humerus can be done also. The following is the illustration of human anatomy planes from Wikipedia.

The elbow has only 1 degree of freedom- flexing and extending. The wrist has 3 degrees of freedom -up and down, and side to side, and rotation of forearm with respective to the axis along ulna. Therefore, the human arm has a total of 7 degrees of freedom.

Byte Stuffing

I occasionally need to write programs to send and receive data bytes from one device to another. That is why I arbitrarily choose a simple variant of byte stuffing methods to build frames to send and receive data. To delimit the frame, control characters - 0x02 and 0x03- are defined as start of text (STX) and end of text (ETX) respectively. For error detection, exclusive-or of data bytes is appended after the ETX as a checksum. If you need better error detection, CRC as described at

CRC Calculation in VB and C

can also be used.

For example, if we want to send two bytes of data -

0x30 0x31,

the resulting frame will be

0x02 0x30 0x31 0x03 0x01

,where 0x02 at the start is added as STX, followed by data bytes and the byte before the last one, 0x03, is added as ETX. Since the exclusive-or of data bytes, 0x02^0x03, is 0x01, it is appended at the end as checksum. How can we send data that contains 0x02 or 0x03 which were already used as control characters? We need to define another control character 0x10 as Data Link Escape (DLE) to mark data that are not control characters. As an another example, let us build a frame for five data bytes -

0x30 0x02 0x65 0x10 0x03.

We will do byte stuffing by putting DLE in front of every data byte that conflicts with STX, ETX, or DLE. And

0x02 0x30 0x10 0x02 0x65 0x10 0x10 0x10 0x03 0x03 0x44

will be the resulting frame. I have developed a few programs in C and LabVIEW. Example programs can be downloaded at the following links.

Byte Stuffing on GitHub


The following is the C++ code to build, send and receive a frame.

// Byte stuffing- sending and receiving frames
// Author: Yan Naing Aye

#ifndef  FRAME_H
#define FRAME_H

#include 
#define STX 0x02
#define ETX 0x03
#define DLE 0x10

#define TX_BUF_SIZE 128
#define RX_BUF_SIZE 128
enum RX_STATE { IGNORE,RECEIVING,ESCAPE,RXCRC1,RXCRC2 };
//-----------------------------------------------------------------------------
class Frame {    
    RX_STATE rState;
protected:
    int TxN;//number of transmitting bytes
    int RxN;//number of receiving bytes
    char tb[TX_BUF_SIZE];//transmit buffer
    char rb[RX_BUF_SIZE];//receiving data
public:
    Frame();
    int setTxFrame(char* d,int n);
    unsigned int CRC16CCITT_Calculate(char* s,unsigned char len,unsigned int crc);
    int getTxN();
    int getRxN();
    int receiveRxFrame(char c);//get receiving frame from received char
    char* getTxBuf();
    char* getRxBuf();
};
//-----------------------------------------------------------------------------
Frame::Frame():TxN(0),RxN(0),rState(IGNORE){}
//-----------------------------------------------------------------------------
char* Frame::getTxBuf(){
    return tb;
}
//-----------------------------------------------------------------------------
char* Frame::getRxBuf(){
    return rb;
}
//-----------------------------------------------------------------------------
//Prepare transmitting frame
int Frame::setTxFrame(char* d,int n)
{
    unsigned int txcrc=0xFFFF;//initialize crc
    char c;
    int i=0,j=0;
    tb[i++]=STX;//start of frame
    for(j=0;j < n;j++) {
        c=d[j];
        if((c==STX)||(c==ETX)||(c==DLE)) tb[i++]=(DLE);
        tb[i++]=c;
    }
    tb[i++]=(ETX);//end of frame

    txcrc=CRC16CCITT_Calculate(d,n,txcrc);//calculate crc
    tb[i++]=txcrc & 0xFF;
    tb[i++]=(txcrc >> 8) & 0xFF;
    TxN=i;
    return TxN;
}
//-----------------------------------------------------------------------------
//Inputs
//s : pointer to input char string
//len: string len (maximum 255)
//crc: initial CRC value

//Output
//Returns calculated CRC
unsigned int Frame::CRC16CCITT_Calculate(char* s,unsigned char len,unsigned int crc)
{
    //CRC Order: 16
    //CCITT(recommendation) : F(x)= x16 + x12 + x5 + 1
    //CRC Poly: 0x1021
    //Operational initial value:  0xFFFF
    //Final xor value: 0
    unsigned char i,j;
    for(i=0;i < len;i++,s++) {
        crc^=((unsigned int)(*s) & 0xFF) << 8;
        for(j=0;j<8;j++) {
            if(crc & 0x8000) crc=(crc << 1)^0x1021;
            else crc <<=1;
        }
    }
    return (crc & 0xFFFF);//truncate last 16 bit
}
//-----------------------------------------------------------------------------
//get number of transmitting bytes
int Frame::getTxN()
{
    return TxN;
}
//-----------------------------------------------------------------------------
//get number of transmitting bytes
int Frame::getRxN()
{
    return RxN;
}
//-----------------------------------------------------------------------------
//process receiving char
int Frame::receiveRxFrame(char c)
{
    static char b;
    unsigned int crc;
    unsigned int rxcrc=0xFFFF;//initialize CRC
    switch(rState){
        case IGNORE:
            if(c==STX) { rState=RECEIVING;RxN=0;}
            break;
        case RECEIVING:
            if(c==STX) { rState=RECEIVING;RxN=0;}
            else if(c==ETX){rState=RXCRC1;}
            else if(c==DLE){ rState=ESCAPE; }
            else { rb[RxN++]=c; }
            break;
        case ESCAPE:
            rb[RxN++]=c; rState=RECEIVING;
            break;
        case RXCRC1:
            b=c; rState=RXCRC2;
            break;
        case RXCRC2:
            rState=IGNORE;
            crc=( (int)c << 8 | ((int)b & 0xFF) ) & 0xFFFF;//get received crc
            rxcrc=CRC16CCITT_Calculate(rb,RxN,rxcrc);//calculate crc
            //printf("crc: %x  rxcrc:%x \n",crc,rxcrc);
            if(rxcrc==crc){return RxN;}//if crc is correct            
            else {RxN=0;}//discard the frame

            break;
    }
    return 0;
}

//-----------------------------------------------------------------------------

//#############################################################################

class Frame2:public Frame {
    char Dt[20];//transmitting data
public:
    Frame2();
    void printTxFrame();
    void printRxFrame();
    void printRxData();
    void setTxData(float x,float y,float z,float b,float t);
};
//-----------------------------------------------------------------------------
Frame2::Frame2():Frame(),Dt(""){}
//-----------------------------------------------------------------------------
//Print out frame content
void Frame2::printTxFrame()
{
    printf("Tx frame buffer: ");
    for(int j=0;j < TxN;j++) printf("%02X ",(unsigned char)tb[j]);
    printf("\n");
}
//-----------------------------------------------------------------------------
//Print out frame content
void Frame2::printRxFrame()
{
    printf("Rx data buffer: ");
    for(int j=0;j < RxN;j++) printf("%02X ",(unsigned char)rb[j]);
    printf("\n");
}
//-----------------------------------------------------------------------------
//Set transmitting data
void Frame2::setTxData(float x,float y,float z,float b,float t)
{
    *(float*)(Dt)=x;
    *(float*)(Dt+4)=y;
    *(float*)(Dt+8)=z;
    *(float*)(Dt+12)=b;
    *(float*)(Dt+16)=t;
    Frame::setTxFrame(Dt,20);
}
//-----------------------------------------------------------------------------
//Print out received data
void Frame2::printRxData()
{
    float x,y,z,b,t;
    x=*(float*)(Dt);
    y=*(float*)(Dt+4);
    z=*(float*)(Dt+8);
    b=*(float*)(Dt+12);
    t=*(float*)(Dt+16);
    printf("Rx data: %f %f %f %f %f \n",x,y,z,b,t);
}
//-----------------------------------------------------------------------------

#endif // FRAME_H

Simple 6 DOF Translation and Rotation Stages

I had the requirement to control a device in 6 degrees of freedom (DOF). There were 3 motorized translation stages to control the position of the tool in 3D space and 3 rotary motorized stages to control its orientation. Although it could be thought as a robot arm having 3 prismatic joints and 3 revolute joints, I did not want to involve complex forward and inverse kinematics for this simple testing equipment.

Here is an example transformation. By arranging rotational axes of rotary motorized stages to intersect at the tool, a simple direct 6 DOF control could be achieved. The following figures show schematic diagram and 3D drawing for the equipment that I designed.

Capacitor selection for a Crystal Oscillator

A typical crystal oscillator circuit and an example layout design are shown in the figures. A crystal manufacturer normally specifies the load capacitance, C_L, needed for a crystal in its datasheet. This value is the load seen by the crystal which is in series with C₁ and C₂, including all parasites. Typically, overall stray capacitance, C_stray, for pin, board, and PCB trace is between 2 pF and 7 pF. The relation between them can be represented by the following formula.


C_L=(C₁ * C₂)/(C₁ + C₂) + C_stray

C₁ and C₂ are normally chosen to have equal value. Therefore,

C₁ = C₂=2 * (C_L - C_stray)

For example, assume that C_stray is 4 pF and crystal load capacitance C_L is 18 pF, C₁ = C₂= 28 pF. The normalized value of 30 pF can be used.

Denavit-Hartenberg Representation of Robots

The D-H model of representation is a simple way of modeling robot links and joints that can be used for any robot configuration. We will have to assign a z-axis and an x-axis for each link. The D-H representation does not use the y-axis at all. Let L_k be the frame associated with link k.

1. Define Z Axis from All Joints

From the axis of joint k+1, define the z^k axis for link k. If the joint is revolute, the z axis is in the direction of rotation as followed by the right hand rule for the rotations. If the joint is prismatic, the z-axis for the joint is along the direction of the linear movement.

2. Define Origins

The intersection of the z^k and z^k-1 axes is selected as the origin of L^k . If they do not intersect, use the intersection of z^k with a common normal between z^k and z^k-1. There is always one line mutually perpendicular to any two skew lines, called common normal, which is the shortest distance between them.

3. Define X-axis

Assign x^k in the direction of the common normal between z^k and z^k-1. If z-axes are intersecting, select x^k to be orthogonal to both z^k and z^k-1 (the direction of the cross-product of the two z-axes). If z^k and z^k-1 are parallel, point x^k away from z^k-1 colinear with the common normal of the previous joint.

4. Define Y-axis

Select y^k to form a right-handed orthonormal coordinate frame L_k.

5. Define the Four Kinematic Parameters

*θ_k is defined as the angle between x^k-1 and x^k axes about the z^k-1 axis. *d_k is the distance between x^k-1 and x^k axes along z^k-1 axis. *a_k is the distance between z^k-1 and z^k axes along x^k axis. *α_k is the angle between z^k-1 and z^k axes about x^k axis.

6. Transferring Frame k-1 to Frame k

*Rotation of L^k-1 about the z^k-1 axis by θ_k will make x^k-1 and x^k parallel to each other. This is true because the common normals a^k-1 and a^k are both perpendicular to z^k-1 axis. *Translation of L^k-1 along z^k-1 axis a distance of d_k will make x^k-1 and x^k colinear. *Translation of L^k-1 along x^k axis a distance of a_k will bring the origins of L^k-1 and L^k together. *Rotation of L^k-1 about x^k axis by α_k will make z^k-1 and z^k axes parallel. At this point, frames L^k-1 and L^k will be exactly the same.



For simulation, a nice free software called Robotassist can be downloaded at www.kinematics.com

Ref:
Introduction to Robotics -Analysis, Control, Applications; Second Edition, Saeed Benjamin Niku
Fundamental of Robotics -Analysis & Control, Robert J. Schilling