MATLAB Examples

Contents

Introduction

In this IoT project, we use a Raspberry Pi, a web cam and ThingSpeak to count cars on a busy highway. We deploy a car-counting algorithm to the Raspberry Pi device, and we analyze and visualize the traffic patterns with ThingSpeak, a data aggregator. This project stores data in channel 38629 on ThingSpeak.

Analytics are everywhere in the Internet of Things (IoT) occurring at:

  1. The edge node
  2. Offline on the desktop
  3. In the cloud

For this project, we show how to develop analytics for the edge device and how to perform exploratory analysis on data collected on the cloud. We also illustrate a simple example of how to perform automated online analysis in the cloud. The example uses ThingSpeak and MATLAB to perform the analyses.

Hardware Setup

We constructed the car counter using a Raspberry Pi 2, and a USB webcam. The webcam, a Logitech HD Pro C920, was mounted on a Monoprice® flexible mini tripod. We placed the camera near a window on the 4th floor of our office building that overlooks Route 9 in Natick, MA. We angled the camera to have a clear view of both sides of the highway. The camera was connected to one USB port of the Raspberry Pi, and we placed a Wi-Pi™ WLAN USB dongle on another USB port. We then connected the Raspberry Pi to the wireless network in the building. The complete parts list is shown below.

Because we did not want to send high-bandwidth video images to the cloud, we chose to detect the vehicles at the edge using the processor on the Raspberry Pi 2. We then send the count value to the data aggregator at an update rate of once every 15 seconds, the maximum data rate allowed by ThingSpeak.

Developing the Car Counting Algorithm

To develop the car counting algorithm, we used Simulink, Image Processing Toolbox, Computer Vision System Toolbox and the Simulink Support Package for Raspberry Pi Hardware. Simulink is a modeling environment that can to automatically generate code that can run on an embedded controller. In this example, Simulink generates code that runs on the Raspberry Pi 2. The Simulink model for the car-counting algorithm is shown below.

To develop the algorithm, we used the external mode capability of Simulink. In this mode, Simulink gathers the video stream from the Raspberry Pi, and the user can view the video on an external monitor using the SDL Video Display block while the algorithm is running, as shown here.

Moving from left to right in the model, we can see the USB webcam connected to the Raspberry Pi captures video with a region of interest selected. Next the vision.ForegroundDetector estimates the foreground pixels of the video sequence captured from the stationary webcam. The vision.ForegroundDetector estimates the background using Gaussian mixture models and produces a foreground mask highlighting foreground objects-in this case, moving cars.

The foreground mask is then post-processed using Median Filter to remove unwanted noise. It is then analyzed using the Blob Analysis block, which computes centroids of the blob containing the cars.

Finally, the Car Counter Block estimates the number of cars traveling in each direction during the 15 second interval. This block breaks up the region of interest into two sub sections along the highway median. It then counts the cars above and below that line. Lastly the data is sent to the data aggregator using the ThingSpeak Write block. In this block, we send a value for eastbound vehicle count and westbound vehicle count. We send both values to ThingSpeak (channel 38629), and store the eastbound value in Field 2 and the westbound value in Field 1.

Accuracy of the Car Counting Algorithm

The count value is an estimate because the algorithm used does not track the cars. As a result, under certain high-volume traffic conditions, such as when traffic is stopped or moving slowly, it may be possible to count some cars twice. To avoid this situation, the algorithm would need to be modified to track cars as they enter and leave the region of interest.

In the following situations, the algorithm might not work as expected:

  1. Because we are not using a night-vision camera, the car-counting algorithm may not give an accurate result at nighttime.
  2. If there is traffic jam on the road segment, cars at or near a standstill become part of background and may not be properly counted.
  3. If there is abrupt change in lighting condition, the algorithm takes couple of seconds to readjust itself.
  4. The algorithm may not count partially occluded cars.

Deploying the Algorithm to the Hardware

Once the algorithm is developed, we can deploy the algorithm to the hardware, and it will run without being connecting to Simulink. This can be accomplished using the Simulink and the Simulink Support Package for Raspberry Pi hardware. For more information on programming the Raspberry Pi with Simulink, see Raspberry Pi Programming Using Simulink (video).

Analyzing Data on ThingSpeak

With the car-counting algorithm deployed on the Raspberry Pi device, we can start to analyze the raw data sent by the Raspberry Pi to ThingSpeak by fetching it from thingspeak.com. To simplify the retrieval of the data from ThingSpeak, we use the functions from the ThingSpeak Support Toolbox, available on MATLAB Central File Exchange.

Reading One Week of Data into MATLAB

We begin by specifying a start date and an end date using a datetime object. Because we are sending data to ThingSpeak once every 15 seconds, we have approximately 40,000 data points to retrieve. ThingSpeak allows only 8000 points in a single read operation, so we create a for loop to gather the data in batches. We then append the traffic and time data from each iteration into two vectors, called alltrafficData and timestamp.

endDate = datetime('8/1/2015', 'InputFormat', 'MM/dd/yyyy');
startDate = datetime('7/25/2015', 'InputFormat', 'MM/dd/yyyy');
% Create date vector
dateVector = startDate: endDate;
% check to see that the last dateVector value is the same as endDate, if
% not append it
if (dateVector(end) ~= endDate)
   dateVector = [dateVector, endDate];
end
alltrafficData = [];
timestamp = [];
% Read data in chunks because ThingSpeak has a limit of 8000 pts per read
for dayCount = 1:length(dateVector)-1
   dateRange = [dateVector(dayCount), dateVector(dayCount+1)];
   [channelData, t] = thingSpeakRead(38629, 'DateRange', dateRange);
   [alltrafficData] = [alltrafficData; channelData ];
   [timestamp] = [timestamp; t];
end

Plotting Raw Traffic Data

Next we plot the data and label the graph. The daily fluctuations in traffic are clearly visible from the raw data. Although harder to see from the raw data, you can observe a different pattern on the weekend days (7/25 and 7/26) compared to the weekdays (7/27-7/31). The weekdays show stronger morning and afternoon peaks due to commuter traffic.

figure
plot(timestamp, alltrafficData)
datetick
xlabel('Date')
ylabel('Vehicle Count per 15 second interval')
grid on
title('Traffic Volume for the week of July 25')
legend('West Bound Traffic','East Bound Traffic')

Looking at Traffic Volume as a Histogram

To get a better idea of the distribution of traffic volume during the week, we plot a histogram of the east bound traffic. In this view, a traffic jam is represented by high vehicle counts and a clear highway is represented by a vehicle count of zero. It is interesting to note that most observations show less than 15 cars passing by on eastbound side of the divided highway which has 3 lanes in each direction. The highest frequency of observations is 0 cars due to the low traffic volume at night.

eastTraffic = alltrafficData(:,2);
westTraffic = alltrafficData(:,1);
figure
histogram(eastTraffic)
grid on
xlim([0 70])
xlabel('East bound vehicles per 15 second interval')
ylabel('Frequency of observations over 1 week')
title('Histogram of traffic flow for the week of July 25')

Looking at Weekly Sums for Eastbound and Westbound Traffic and Plotting

To visualize the volume of traffic each day, we can easily sum up the raw traffic data and use a bar chart to visualize it. Interestingly, we see that the eastbound traffic on this particular weekend was even higher than the weekday traffic.

for i=1:6
eastTraffic(i) = [];
westTraffic(i) = [];
end
% Divide data into 7 chunks (approximates a 24 hour day)
dailysumeast = sum(reshape(eastTraffic, floor(length(alltrafficData)/7),[]));
dailysumwest = sum(reshape(westTraffic, floor(length(alltrafficData)/7),[]));
dates = dateVector; % convert to serial date for bar plot
dates(8) = [];
dates = datenum(dates);
figure
bar(dates,dailysumeast)
grid on
xlabel('Date')
ylabel('Vehicle Count per day')
title('East Bound Traffic Volume for the week of July 25')
datetick
figure
bar(dates, dailysumwest)
grid on
xlabel('Date')
ylabel('Vehicle Count per day')
title('West Bound Traffic Volume for the week of July 25')
datetick

Taking a Deeper Look at Individual Days during the Week and on the Weekend

To gain more insight into our traffic data, we look at traffic on the individual days. For each day, we specify a start time and a stop time.

clear all
startTime{1} = 'July 27, 2015 00:00:00';
stopTime{1} = 'July 27, 2015 23:59:59';
startTime{2} = 'July 28, 2015 00:00:00';
stopTime{2}= 'July 28, 2015 23:59:59';
startTime{3} = 'July 29, 2015 00:00:00';
stopTime{3} = 'July 29, 2015 23:59:59';
startTime{4} = 'July 30, 2015 00:00:00';
stopTime{4} = 'July 30, 2015 23:59:59';
startTime{5} = 'July 31, 2015 00:00:00';
stopTime{5} = 'July 31, 2015 23:59:59';
startTime{6 }= 'July 25, 2015 00:00:00'; % weekend day
stopTime{6} = 'July 25, 2015 23:59:59';
startTime{7}= 'July 26, 2015 00:00:00'; % weekend day
stopTime{7} = 'July 26, 2015 23:59:59';

Fetch Each Individual Day in a Loop

For each of the 7 days, we retrieve the data from ThingSpeak, and we downsample the data to remove short-term fluctuations. We then use the findpeaks function to plot and find the times where the traffic volume is highest, and we make some observations for each day. For simplicity, we look at the eastbound data only.

for jj = 1:7
    % This loop executes once for each day of the week
startDate = datetime(startTime{jj}, 'InputFormat', 'MMMM d, yyyy HH:mm:ss ');
endDate = datetime(stopTime{jj}, 'InputFormat', 'MMMM d, yyyy HH:mm:ss ');
datevector = [startDate, endDate];
[Daily, t] = thingSpeakRead(38629, 'DateRange', datevector);
% Look at east bound traffic only
DailyEast = Daily(:, 2);
timestamp = datetime(t,'ConvertFrom','datenum');
dateAnalyzed = startTime{jj};
dateAnalyzed = {dateAnalyzed(1:(end-8))};

Check Time Difference Between Samples to Verify How Much Data is Missing

We expect data to be sent from the device once every 15 seconds. In reality, the mean distance between samples is only approximately 15 seconds and the maximum distance between samples is just over 30 seconds. This indicates that there are occasional dropped samples. It also explains why the number of samples in a day is less than 24*60*4=5760. Because the sampling is not exactly uniform, we need to choose a way to bin the data when we downsample.

timediff = diff(timestamp);
if jj == 1 % lets look at the first day
dateAnalyzed
maxtimediff = max(timediff)
meantimemdiff = mean(timediff)
end
dateAnalyzed = 

    'July 27, 2015 '


maxtimediff = 

   00:00:29


meantimemdiff = 

   00:00:15

Downsample into 48 Bins of Approximately 30 Minute Chunks of Data and Find Peaks

The raw traffic data is very spiky and hard to visualize. If we want to see what time of day has the highest volume of traffic, we need to look at the data on a time scale larger than 15 seconds. To do this we divide the 24 hour day into 30 minute segments. Each segment begins near the top of the hour and ends at approximately 30 minutes later.

downsamplesize = floor(length(DailyEast)/48);
teastper30 = downsample(t, downsamplesize);
DailyEastper30(1:48) = 0; % pre-allocate
for k = 1:48
DailyEastper30(k) = sum(DailyEast(1+downsamplesize*(k-1):downsamplesize*k));
end
teastper30(1) = []; % start first bin at 12:30 am
timestampPer30=datetime(teastper30,'ConvertFrom','datenum');
% Find peaks and their times (locations)
[peaks,location] = findpeaks(DailyEastper30,  'Threshold',100, 'MinPeakHeight', 1100);
% Plot peaks
figure
findpeaks(DailyEastper30, datenum(timestampPer30),'Threshold',100, 'MinPeakHeight', 1100)
datetick
xlabel('Time of Day')
ylabel('Vehicle count per 30 minutes')
title(strcat('Peak volume on ', {' '}, dateAnalyzed))
dateAnalyzed
peaktimes = timestampPer30(location)
DailyVolume = sum(DailyEast)
dateAnalyzed = 

    'July 27, 2015 '


peaktimes = 

   27-Jul-2015 06:25:26
   27-Jul-2015 07:25:22
   27-Jul-2015 08:25:19
   27-Jul-2015 16:25:06
   27-Jul-2015 18:25:18
   27-Jul-2015 19:24:43
   27-Jul-2015 20:24:04


DailyVolume =

       43667

dateAnalyzed = 

    'July 28, 2015 '


peaktimes = 

   28-Jul-2015 06:57:01
   28-Jul-2015 14:26:27
   28-Jul-2015 15:58:15
   28-Jul-2015 16:59:26
   28-Jul-2015 20:57:43


DailyVolume =

       50327

dateAnalyzed = 

    'July 29, 2015 '


peaktimes = 

   29-Jul-2015 06:28:17
   29-Jul-2015 07:28:24
   29-Jul-2015 19:25:46
   29-Jul-2015 20:55:09


DailyVolume =

       37850

dateAnalyzed = 

    'July 30, 2015 '


peaktimes = 

   30-Jul-2015 06:26:27
   30-Jul-2015 07:56:27
   30-Jul-2015 14:26:31
   30-Jul-2015 16:26:16
   30-Jul-2015 18:57:52
   30-Jul-2015 20:26:38


DailyVolume =

       42146

dateAnalyzed = 

    'July 31, 2015 '


peaktimes = 

   31-Jul-2015 05:27:18
   31-Jul-2015 06:26:33
   31-Jul-2015 16:58:57
   31-Jul-2015 19:29:46
   31-Jul-2015 20:29:22


DailyVolume =

       45741

dateAnalyzed = 

    'July 25, 2015 '


peaktimes = 

   25-Jul-2015 08:55:31
   25-Jul-2015 09:54:42
   25-Jul-2015 10:53:56
   25-Jul-2015 14:22:14
   25-Jul-2015 15:22:22
   25-Jul-2015 16:21:43
   25-Jul-2015 17:20:50
   25-Jul-2015 19:51:11


DailyVolume =

       68923

dateAnalyzed = 

    'July 26, 2015 '


peaktimes = 

   26-Jul-2015 10:25:52
   26-Jul-2015 12:24:50
   26-Jul-2015 14:53:03
   26-Jul-2015 16:51:46
   26-Jul-2015 19:49:26


DailyVolume =

       60434

Observations for Monday, July 27

We can clearly see the morning and evening peak-traffic times for this particular Monday morning. There are three 30 minute windows with peaks in the morning rush and 4 peaks in the afternoon rush. On this particular day, there was a small peak at lunch time, but it was lower than the threshold we set for the peak finder.

Observations for Tuesday, July 28

Again we can clearly see the morning and evening peak-traffic times for this Tuesday morning. There is 1 large peak in the morning and 5 peaks in the afternoon. One of the afternoon peaks is before 3 p.m. and may be the result of an accident or construction.

Observations for Wednesday, July 29

On Wednesday, we can clearly see the morning and evening peak-traffic conditions although unlike on Monday we see only 2 peaks in the morning in two peaks in the evening on this particular day. Note that the evening peak-travel times peak at a higher volume (more than 3000 vehicles in 30 minutes) than the morning peak-travel volumes.

Observations for Thursday, July 30

The observed pattern on Thursday looks similar to Wednesday although we see two extra peaks in the early afternoon. This could be due to an accident on the roadway.

Observations for Friday, July 31

Finally on Friday we see 2 peaks in the morning rush and 3 peaks in the afternoon rush. The largest peak is approximately 5000 vehicles per 30 minutes and occurs between 7 and 7:30 p.m. This peak may be higher than other days due to people travelling for a summer weekend trip.

Weekend observations for Saturday, July 25

On Saturday, we see a totally different traffic pattern with many peaks scattered almost evenly throughout the day. Although there are no large rush-hour peaks, the average volume is quite high from 9 a.m. through 8 p.m., with many time periods seeing more than 2000 vehicles in a half hour. These results are not surprising because this particular highway provides access to multiple malls and car dealerships.

Weekend observations for Sunday, July 26

On Sunday, we see yet a different pattern with some heavy traffic in mid-morning at around 9:55 a.m. to 10:25 a.m. and then another peak in the evening at between 7:20 and 7:50 p.m.

end

Online Analysis: Calculating and Scheduling MATLAB Code inside ThingSpeak

So far, we’ve developed a car-counting algorithm that was deployed onto the Raspberry Pi 2 hardware, and we've programmed it to send its results to the ThingSpeak data aggregation service. Next we, brought the data back into MATLAB to do some offline analysis of the daily and weekly traffic patterns. We used the ThingSpeak Support Toolbox, available on MATLAB Central File Exchange to retrieve the data for analysis.

But what if we want to compute the daily volume of traffic and send that data to a channel on ThingSpeak? It’s very simple to do the calculation offline, but ThingSpeak with its MATLAB integration enables us to do this online too. To calculate the daily volume of traffic, we use the MATLAB Analysis App available on Thingspeak.com. We write code to fetch the previous day’s raw traffic data and we sum up the results. We then use the thingSpeakWrite command to send the data to a new channel called Daily Volume.

Finally, we schedule our analysis to run once a day just after midnight using the ThingSpeak TimeControl App.

Viewing the Calculated Data Online

We have now created a channel that reports the daily volume of cars counted. Because this channel is on ThingSpeak and the MATLAB Analysis App is built in, we can look at this analysis at any time using any web browser on a PC or mobile device. To see the calculated value, go to channel 51671 on ThingSpeak.

Conclusions and Acknowledgements

Analytics are everywhere in IoT! In this article, we demonstrated how to develop a car-counting algorithm that can be deployed onto Raspberry Pi edge node and sends data to ThingSpeak. We showed how to use MATLAB for offline analysis by retrieving data from ThingSpeak and analyzing and visualizing daily and weekly traffic patterns. And we showed how to perform online analysis inside the ThingSpeak web service by using the MATLAB Analysis App to create a daily traffic-volume channel that automatically updates once a day.

I would like to thank Mohammad Khan for his contribution to this project and this article.

Copyright 2015, The MathWorks, Inc.