A Recursive Least Square harmonic estimation algorithm is implemented in Arduino Due microcontroller. A Laboratory setup to generate a test signal, i.e. a modified sine wave inverter was used to generate a 24 volt AC signal, which is sensed by the Microcontroller through a signal conditioning circuit and then the microcontroller finds out the harmonic components and communicates them to the computer which are plotted in MATLAB.

To make the AC Voltage Signal Get Accepted by the ADC of the Microcontroller we need the following circuits.

Voltage Supply Circuit (230 Volt AC to 5 Volt DC):

  1. 230 Volt to 6 Volt Step-Down Transformer: It steps down the voltage to 6 Volt to be used by our Power Supply Circuit for the Microcontroller as well as the Voltage Offset Circuit.
  2. Diode Bridge Rectifier: A bridge Rectifier made of Four Diodes rectify the 6 Volt AC of the Transformer secondary, i.e. Only Positive Parts of the signals i.e. DC.
  3. Smoothing Capacitor: A 1000uF 25 Volt Smoothing capacitor smooths out the Voltage coming from the Bridge Rectifier.
  4. CD7805 Voltage Regulator: A Voltage regulator eliminates any ripple present in the Voltage supply after the capacitor so that a maximum allowable ripple of .005 Volt is present, i.e. the Output is close to Pure DC suitable for the sensitive Microcontroller.

Measurement Signal Pre-Processing Circuit:

  1. The 9 Volt Transformer: Since we will be measuring 230 Volt AC it has to be stepped down to a suitable voltage range to be operated on and also the Transformer Provides Sufficient Isolation necessary to keep ourselves as well as the sophisticated electronics from harm’s way.
  2. Voltage Divider (1/10): The output of the Transformer, when unloaded can give up to 14-17 volt Peak. Hence it has to be further reduced, for this purpose we have used a Voltage divider which divides it 10 times which is the ratio of the resistances 10K and 100K.
  3. Voltage Offset Adding Circuit: The voltage after the voltage divider has both a positive and negative peak, which is not suitable for our ADC. Therefore, we add an OFFSET using the Rail to Rail swing OPAMP LM358. The OPAMP functions as a voltage follower, and buffers the mid-point voltage that appears at the junction of R3 & R4. This significantly reduces the impedance of the voltage source, resulting in enhanced performance.
  4. The Final Voltage Divider: Since the original circuit gives out a signal in the range of 5 Volt, but since we have used a microcontroller, which has an ADC capable of only measuring a signal in the range of 3.3 Volt, a final Voltage divider has been used to keep the voltage in range suitable for the ADC.

 

The Experimental Setup looks Like the Following:

 

The Output waveform of the PWM inverter whose Harmonics are to be estimated

Estimated Amplitudes Of The Harmonics

Estimated Phases Of The Harmonics

The MATLAB Code for SImulation:

Matlab M Code:
  1. clear all;
  2. clc;
  3. %tic
  4. dt=.001;
  5. f=50;
  6. w=2*pi*f;
  7. ph1=4*pi/9;
  8. ph3=pi/3;
  9. ph5=pi/4;
  10. ph7=pi/5;
  11. ph11=pi/6;
  12. adc=0.5;
  13. alpha=5;
  14. a1=1.5;a3=0.5;a5=0.2;a7=0.15;a11=0.1;
  15. for i=1:100
  16. t(i)=dt*i;
  17. y(i)=a1*sin(w*dt*i+ph1)+a3*sin(3*w*dt*i+ph3)+a5*sin(5*w*dt*i+ph5)+a7*sin(7*w*dt*i+ph7)+a11*sin(11*w*dt*i+ph11)+adc*exp(-5*dt*i)+.01*randn(1,1);
  18.  
  19. h11(i)=[sin(w*dt*i)];
  20. h12(i)=[cos(w*dt*i)];
  21. h21(i)=[sin(3*w*dt*i)];
  22. h22(i)=[cos(3*w*dt*i)];
  23. h31(i)=[sin(5*w*dt*i)];
  24. h32(i)=[cos(5*w*dt*i)];
  25. h41(i)=[sin(7*w*dt*i)];
  26. h42(i)=[cos(7*w*dt*i)];
  27. h51(i)=[sin(11*w*dt*i)];
  28. h52(i)=[cos(11*w*dt*i)];
  29. h61(i)=i;
  30. h62(i)=-dt*i;
  31. end
  32. R=[h11 ;h12 ; h21 ;h22; h31; h32; h41; h42; h51; h52; h61; h62]';
  33. P=diag(100:111);
  34.  
  35. for k=1:length(y)
  36. sitrans(k,:)=R(k,:);% dimenson of 100*6
  37. si=sitrans';% dimenson of 6*100
  38. end
  39. th=[0;0;0;0;0;0;0;0;0;0;0;0];
  40. for m=1:length(y)
  41. P=P-[(P*si(:,m)*sitrans(m,:)*P)/(1+sitrans(m,:)*P*si(:,m))];
  42. %P=P-[(P*R'*R*P)/(1+R'*R*P)];
  43. k=P*si(:,m);
  44. %e(m)=y(m)-sitrans(m,:)*th;
  45. e(m)=y(m)-R(m,:)*th;
  46. th=th+k*e(m);
  47.  
  48. yo1(m)=[h11(m) h12(m)]*[th(1,1); th(2,1)];
  49. yo3(m)=[h21(m) h22(m)]*[th(3,1); th(4,1)];
  50. yo5(m)=[h31(m) h32(m)]*[th(5,1); th(6,1)];
  51. yo7(m)=[h41(m) h42(m)]*[th(7,1); th(8,1)];
  52. yo11(m)=[h51(m) h52(m)]*[th(9,1); th(10,1)];
  53. yo(m)=sitrans(m,:)*th;
  54.  
  55. %%rmse(m)=sqrt(0.5*(e(m)).^2);
  56. %%mse(m)=(0.01*(e(m).^2));
  57.  
  58. amp1(m)=sqrt((th(1,1))^2+(th(2,1))^2);
  59. amp3(m)=sqrt((th(3,1))^2+(th(4,1))^2);
  60. amp5(m)=sqrt((th(5,1))^2+(th(6,1))^2);
  61. amp7(m)=sqrt((th(7,1))^2+(th(8,1))^2);
  62. amp11(m)=sqrt((th(9,1))^2+(th(10,1))^2);
  63.  
  64. phased1(m)=atand(th(2,1)/(th(1,1)));
  65. phased3(m)=atand(th(4,1)/(th(3,1)));
  66. phased5(m)=atand(th(6,1)/(th(5,1)));
  67. phased7(m)=atand(th(8,1)/(th(7,1)));
  68. phased11(m)=atand(th(10,1)/(th(9,1)));
  69.  
  70. end
  71. % toc
  72.  

The Arduino Microcontroller Code for Realtime Estimation



Arduino Code:
  1. #include "Eigen312.h" // Calls main Eigen3.1.2 matrix class library
  2. #include "LU" // Calls inverse, determinant, LU decomp., etc.
  3. using namespace Eigen; // Eigen related statement; simplifies syntax for declaration of matrices
  4.  
  5. void print_mtxf(const Eigen::MatrixXf& K);
  6.  
  7. float dt = .001;
  8. int f = 50;
  9. float pi = 3.143;
  10. float w = 2 * pi * f;
  11. float ph1 = 4 * pi / 9, ph3 = pi / 3, ph5 = pi / 4, ph7 = pi / 5, ph11 = pi / 6;
  12. float adc = 0.5;
  13. int alpha = 5;
  14. float a1 = 1.5, a3 = 0.5, a5 = 0.2, a7 = 0.15, a11 = 0.1;
  15. int i;
  16. int j;
  17. int m;
  18. MatrixXf y(1, 100); // Note: without "using namespace Eigen", declaration would be: Eigen::MatrixXf H(6,6);
  19. //MatrixXf volt(1,100);
  20. MatrixXf actual_data(1, 100);
  21. void setup() {
  22.  
  23. Serial.begin(115200);
  24.  
  25. // DECLARE MATRICES
  26. //--------------------
  27. MatrixXf t(1, 100); // Produces 6x6 float matrix class
  28. //MatrixXf y(1,100); // Note: without "using namespace Eigen", declaration would be: Eigen::MatrixXf H(6,6);
  29.  
  30. MatrixXf h11(1, 100);
  31. MatrixXf h12(1, 100);
  32. MatrixXf h21(1, 100);
  33. MatrixXf h22(1, 100);
  34. MatrixXf h31(1, 100);
  35. MatrixXf h32(1, 100);
  36. MatrixXf h41(1, 100);
  37. MatrixXf h42(1, 100);
  38. MatrixXf h51(1, 100);
  39. MatrixXf h52(1, 100);
  40.  
  41. MatrixXf h61(1, 100);
  42. MatrixXf h62(1, 100);
  43.  
  44. MatrixXf Rt(12, 100);
  45. MatrixXf R(100, 12);
  46.  
  47. MatrixXf P(12, 12);
  48.  
  49. MatrixXf sitrans(100, 12);
  50. MatrixXf si(12, 100);
  51.  
  52. MatrixXf th(12, 1);
  53. MatrixXf k(12, 1);
  54. MatrixXf e(1, 100);
  55.  
  56. MatrixXf yo1(1, 100);
  57. MatrixXf yo3(1, 100);
  58. MatrixXf yo5(1, 100);
  59. MatrixXf yo7(1, 100);
  60. MatrixXf yo11(1, 100);
  61.  
  62. MatrixXf yo(1, 100);
  63.  
  64. MatrixXf amp1(1, 100);
  65. MatrixXf amp3(1, 100);
  66. MatrixXf amp5(1, 100);
  67. MatrixXf amp7(1, 100);
  68. MatrixXf amp11(1, 100);
  69.  
  70. MatrixXf phased1(1, 100);
  71. MatrixXf phased3(1, 100);
  72. MatrixXf phased5(1, 100);
  73. MatrixXf phased7(1, 100);
  74. MatrixXf phased11(1, 100);
  75.  
  76. // Print Result
  77. //----------------------------
  78. //print_mtxf(t); // Print Matrix Result (passed by reference)
  79. get_voltage_samples();
  80.  
  81. for (int x = 1; i < 101; i++)
  82. {
  83. t(0, i - 1) = dt * i;
  84.  
  85. //y(0, i - 1) = a1 * sin(w * dt * i + ph1) + a3 * sin(3 * w * dt * i + ph3) + a5 * sin(5 * w * dt * i + ph5) + a7 * sin(7 * w * dt * i + ph7) + a11 * sin(11 * w * dt * i + ph11) + adc * exp(-5 * dt * i) + .01 * random(-1, 1) * i;
  86.  
  87.  
  88.  
  89. h11(0, i - 1) = sin(w * dt * i);
  90. h12(0, i - 1) = cos(w * dt * i);
  91. h21(0, i - 1) = sin(3 * w * dt * i);
  92. h22(0, i - 1) = cos(3 * w * dt * i);
  93. h31(0, i - 1) = sin(5 * w * dt * i);
  94. h32(0, i - 1) = cos(5 * w * dt * i);
  95. h41(0, i - 1) = sin(7 * w * dt * i);
  96. h42(0, i - 1) = cos(7 * w * dt * i);
  97. h51(0, i - 1) = sin(11 * w * dt * i);
  98. h52(0, i - 1) = cos(11 * w * dt * i);
  99.  
  100. h61(0, i - 1) = i;
  101. h62(0, i - 1) = -dt * i;
  102.  
  103. }
  104.  
  105. Serial.println("t Matrix is");
  106. print_mtxf(t);
  107.  
  108. Serial.println("Actual Data Read from ADC is");
  109. print_mtxf(actual_data);
  110.  
  111. //Serial.println("Mapped Voltage Matrix is");
  112. //print_mtxf(volt);
  113.  
  114. Serial.println("y Matrix is");
  115. print_mtxf(y);
  116.  
  117. Rt.row(0) = h11;
  118. Rt.row(1) = h12;
  119. Rt.row(2) = h21;
  120. Rt.row(3) = h22;
  121. Rt.row(4) = h31;
  122. Rt.row(5) = h32;
  123. Rt.row(6) = h41;
  124. Rt.row(7) = h42;
  125. Rt.row(8) = h51;
  126. Rt.row(9) = h52;
  127. Rt.row(10) = h61;
  128. Rt.row(11) = h62;
  129.  
  130.  
  131. Serial.println("Rt Matrix is");
  132. print_mtxf(Rt);
  133. R = Rt.transpose();
  134. Serial.println("R Matrix is");
  135. print_mtxf(R);
  136.  
  137. P.setZero();
  138. for (i = 0; i < 12; i++) {
  139. for (j = 0; j < 12; j++) {
  140. if (i == j) {
  141. P(i, j) = 100 + i;
  142. }
  143. }
  144. }
  145. Serial.println("P matrix");
  146. print_mtxf(P);
  147.  
  148. sitrans = R;
  149. si = sitrans.transpose();
  150.  
  151. Serial.println("si matrix is");
  152. print_mtxf(si);
  153.  
  154. th.setZero();
  155.  
  156. for (m = 0; m < 100; m++) {
  157. //P=P-[(P*si(:,m)*sitrans(m,:)*P)/(1+sitrans(m,:)*P*si(:,m))];
  158. // temp1 temp2
  159. P = P - ((P * si.col(m) * sitrans.row(m) * P) / (1 + sitrans.row(m) * P * si.col(m)));
  160. k = P * si.col(m);
  161. //e(m)=y(m)-R(m,:)*th;
  162. e(0, m) = y(0, m) - (R.row(m) * th)(0, 0);
  163. // temp1=R.row(m)*th;
  164. // e(0,m)= y(0,m) - temp1(0,0);
  165. th = th + k * e(0, m);
  166.  
  167. //yo1(m)=[h11(m) h12(m)]*[th(1,1); th(2,1)];
  168. MatrixXf temp1(1, 2);
  169. MatrixXf temp2(2, 1);
  170.  
  171. temp1(0, 0) = h11(0, m);
  172. temp1(0, 1) = h12(0, m);
  173.  
  174. temp2(0, 0) = th(0, 0);
  175. temp2(1, 0) = th(1, 0);
  176.  
  177. yo1(0, m) = (temp1 * temp2)(0, 0);
  178.  
  179. //yo3(m)=[h21(m) h22(m)]*[th(3,1); th(4,1)];
  180. temp1(0, 0) = h21(0, m);
  181. temp1(0, 1) = h22(0, m);
  182.  
  183. temp2(0, 0) = th(2, 0);
  184. temp2(1, 0) = th(3, 0);
  185.  
  186. yo3(0, m) = (temp1 * temp2)(0, 0);
  187.  
  188. //yo5(m)=[h31(m) h32(m)]*[th(5,1); th(6,1)];
  189. temp1(0, 0) = h31(0, m);
  190. temp1(0, 1) = h32(0, m);
  191.  
  192. temp2(0, 0) = th(4, 0);
  193. temp2(1, 0) = th(5, 0);
  194.  
  195. yo5(0, m) = (temp1 * temp2)(0, 0);
  196.  
  197. //yo7(m)=[h41(m) h42(m)]*[th(7,1); th(8,1)];
  198. temp1(0, 0) = h41(0, m);
  199. temp1(0, 1) = h42(0, m);
  200.  
  201. temp2(0, 0) = th(6, 0);
  202. temp2(1, 0) = th(7, 0);
  203.  
  204. yo7(0, m) = (temp1 * temp2)(0, 0);
  205.  
  206. //yo11(m)=[h51(m) h52(m)]*[th(9,1); th(10,1)];
  207. temp1(0, 0) = h51(0, m);
  208. temp1(0, 1) = h52(0, m);
  209.  
  210. temp2(0, 0) = th(8, 0);
  211. temp2(1, 0) = th(9, 0);
  212.  
  213. yo11(0, m) = (temp1 * temp2)(0, 0);
  214.  
  215.  
  216. yo(0, m) = (sitrans.row(m) * th)(0, 0);
  217.  
  218. //amp1(m)=sqrt((th(1,1))^2+(th(2,1))^2);
  219. //amp1(0,m)=(th(0,0).square()+th(1,0).square()).sqrt();
  220. amp1(0, m) = sqrt(th(0, 0) * th(0, 0) + th(1, 0) * th(1, 0));
  221. amp3(0, m) = sqrt(th(2, 0) * th(2, 0) + th(3, 0) * th(3, 0));
  222. amp5(0, m) = sqrt(th(4, 0) * th(4, 0) + th(5, 0) * th(5, 0));
  223. amp7(0, m) = sqrt(th(6, 0) * th(6, 0) + th(7, 0) * th(7, 0));
  224. amp11(0, m) = sqrt(th(8, 0) * th(8, 0) + th(9, 0) * th(9, 0));
  225.  
  226. //phased1(m)=atand(th(2,1)/(th(1,1)));
  227. phased1(0, m) = atan(th(1, 0) / (th(0, 0)));
  228. phased3(0, m) = atan(th(3, 0) / (th(2, 0)));
  229. phased5(0, m) = atan(th(5, 0) / (th(4, 0)));
  230. phased7(0, m) = atan(th(7, 0) / (th(6, 0)));
  231. phased11(0, m) = atan(th(9, 0) / (th(8, 0)));
  232.  
  233. }
  234. //Convert raian phase matrices into degree
  235. phased1 = phased1 * (180 / pi);
  236. phased3 = phased3 * (180 / pi);
  237. phased5 = phased5 * (180 / pi);
  238. phased7 = phased7 * (180 / pi);
  239. phased11 = phased11 * (180 / pi);
  240.  
  241. Serial.println("yo1 Matrix is");
  242. print_mtxf(yo1);
  243. Serial.println("yo3 Matrix is");
  244. print_mtxf(yo3);
  245. Serial.println("yo5 Matrix is");
  246. print_mtxf(yo5);
  247. Serial.println("yo7 Matrix is");
  248. print_mtxf(yo7);
  249. Serial.println("yo11 Matrix is");
  250. print_mtxf(yo11);
  251.  
  252. Serial.println("amp1 Matrix is");
  253. print_mtxf(amp1);
  254. Serial.println("amp3 Matrix is");
  255. print_mtxf(amp3);
  256. Serial.println("amp5 Matrix is");
  257. print_mtxf(amp5);
  258. Serial.println("amp7 Matrix is");
  259. print_mtxf(amp7);
  260. Serial.println("amp11 Matrix is");
  261. print_mtxf(amp11);
  262.  
  263.  
  264. Serial.println("phased1 Matrix is");
  265. print_mtxf(phased1);
  266. Serial.println("phased3 Matrix is");
  267. print_mtxf(phased3);
  268. Serial.println("phased5 Matrix is");
  269. print_mtxf(phased5);
  270. Serial.println("phased7 Matrix is");
  271. print_mtxf(phased7);
  272. Serial.println("phased11 Matrix is");
  273. print_mtxf(phased11);
  274. }
  275.  
  276.  
  277. void loop() {
  278. // put your main code here, to run repeatedly:
  279.  
  280. }
  281.  
  282. void get_voltage_samples()
  283. {
  284. for (int i = 1; i < 101; i++)
  285. {
  286. float sample = analogRead(A0);
  287. actual_data(0, i - 1) = sample;
  288. float sample_mapped_230 = map(sample, 108, 426, -320.00, 320.00);
  289. //volt(0, i - 1)=sample_mapped_230;
  290. //float sample_mapped_ex=map_double(sample,108,426,-2.6000,2.6000);
  291. y(0, i - 1) = sample_mapped_230;
  292. delay(1);
  293.  
  294. //delayMicroseconds(100);
  295. }
  296. }
  297.  
  298. float map_double(double x, double in_min, double in_max, double out_min, double out_max)
  299. {
  300. return (x - in_min) * (out_max - out_min) / (in_max - in_min) + out_min;
  301. }
  302.  
  303.  
  304. // PRINT MATRIX (float type)
  305. //-----------------------------
  306. void print_mtxf(const Eigen::MatrixXf& X)
  307. {
  308. int i, j, nrow, ncol;
  309.  
  310. nrow = X.rows();
  311. ncol = X.cols();
  312.  
  313. Serial.print("nrow: "); Serial.println(nrow);
  314. Serial.print("ncol: "); Serial.println(ncol);
  315. Serial.println();
  316.  
  317. for (i = 0; i < nrow; i++)
  318. {
  319. for (j = 0; j < ncol; j++)
  320. {
  321. Serial.print(X(i, j), 6); // print 6 decimal places
  322. Serial.print(", ");
  323. }
  324. Serial.println();
  325. }
  326. Serial.println();
  327. }