Differences

This shows you the differences between two versions of the page.

--- wind_sensor:wind_sensor_documentation [2015/11/25 21:05]
jeremygg [Preprocessing the Data]
+++ wind_sensor:wind_sensor_documentation [2021/09/19 21:59] (current)
@@ Line 15: / Line 15: @@
    - squaring a function in time {x(t) -> x(t)<sup>2</sup>} is convolution in frequency domain {X(f) -> X(f)∗X(f)}. Also, sound is a real base band signal, so there for it is a low passed signal which is symmetrical around 0. (Squaring is similar to taking the absolute value. Either method could work)
    - A convolution between the same signal will be symmetrical and have it's peak in the middle. In this case it is at 0, and X(0) is the DC part which is the total energy in the signal.
-   - You can find the frequency scaling depending on your sampling frequency. If sampling frequency is ω<sub>s</sub>, the frequency vector for the fft would go from -ω<sub>s</sub>/2 to ω<sub>s</sub>.
+   - You can find the frequency scaling depending on your sampling frequency. If sampling frequency is ω<sub>s</sub>, the frequency vector for the fft would go from -ω<sub>s</sub>/2 to ω<sub>s</sub>/2.
 {{:wind_sensor:untitled2.png?nolink&400|}}
@@ Line 42: / Line 42: @@
    - We noted that the output graph looked like a square root function, and after squaring, it looked like a linear line,so we predicted on the square of the output, and then we will square root the prediction.
-(Insert Code Snippet)
+    % For actual fitting you could use either the moving average or regular average (CT or RT).
+    % We noticed it look like a trig function, so we just used sin and cos to
+    % fit it.
+    %Since we our going from output to angle, instead of angle to output, we
+    %can't make an explicit function since trig function is not aninvertible
+    %function between 0 and 2pi. This is why we are going to overlap from 0 to
+    %pi/2
+    theta = mod(theta,360);                        % Degrees are circular, so deg = 360*k + deg
+    phi = mod(theta,180);
+    phi = (theta>180).*(180-phi)+phi.*(theta<180)+180.*(theta==180);
+    figure(4)
+    hold on
+    scatter(CT,phi,'g')
+    plot(CT,pred)
+    Label('Compare prediction to actual','peak detector','theta');
+    legend('Measured','Predicted')
+{{:wind_sensor:angpred.jpg|}}
-(Insert Linear Model for Angles, Wind Sensor)
+    %We note we are off at low values, so we decided to do a piecewise
+    %prediction, 2 different models between first 4 samples and last samples.
+    %Note the 4rth sample happens at output < 13.
+    bp = 4;                                                   % Break point is at fourth sample
+    model = polyfit(CT(1:bp),Speed(1:bp).^2,1);               % fit a model up to Nth degree till break point
+    pred = sqrt(polyval(model,CT(1:bp)));                     % evaluate model with input, then squared
+    model = polyfit(CT(bp+1:end),Speed(bp+1:end).^2,1);       % fit a model up to Nth degree after break point
+    pred = [pred sqrt(polyval(model,CT(bp+1:end)))];          % evaluate model with input, then squared
+    figure(5)
+    plot(CT,Speed,CT,pred)
+    Label('Compare prediction to actual','peak detector','mph');
+    legend('Measured','Predicted');
+{{:wind_sensor:amppred2.jpg|}}
 Note we tampered but not fully explored in the frequency domain. Looking at these angles, or amplitudes in general, as speed gets higher the tone gets higher which means the frequency gets higher. (You could use the Find_peak function needed in EE415)This is shown with the impulses getting shifted outward as amplitudes gets higher.
-(close up picture of impulses)
+{{:wind_sensor:imp.jpg|}}
 This symmetrical impulses at frequency -f<sub>c</sub> and f<sub>c</sub>, corresponds to the Fourier transform of a cosine. Notice how majority of the weight of the signal has most of its energy and magnitude at this impulse. If you also look closely, it is technically also added with a band passed signal which could be due to the microphone specifications. The impulse could be the sound of the fan, but speculation from using different fans and wind speeds lead us to believe that this impulse is the fundamental frequency of the noise of the wind. If you are low passing and sampling, try to stay above the cut off frequencies and stay way above tolerance of that impulse.
@@ Line 60: / Line 94: @@
 What is needed is a memory-less way to update the mean of the current signal.
-In this case we do a moving average type of equation. Assume you have a signal of N samples with samples k<sub>1</sub> to k<sub>N</sub> with a mean of μ. If you would only like to take the average of the past N samples, after adding sample k<sub>N+1</sub>, the average would be (k<sub>N+1</sub>+Nμ-k<sub>1</sub>)/N. If the variance of the k's are low and/or the N is high enough, you can treat k<sub>1</sub> practically as μ. You can calculate for the new mean, lets call it μ[n+1] = (k<sub>n</sub>+Nμ[n]-μ[n])/N = μ[n+1] (k<sub>n</sub>)/N+((N-1)μ[n])/N If the signal ends up changing  from its original mean to a different value, the steady state of this method will converge to the newer mean. Our implementation is below:
+In this case we do a moving average type of equation. Assume you have a signal of N samples with samples k<sub>1</sub> to k<sub>N</sub> with a mean of μ. If you would only like to take the average of the past N samples, after adding sample k<sub>N+1</sub>, the average would be (k<sub>N+1</sub>+Nμ-k<sub>1</sub>)/N. If the variance of the k's are low and/or the N is high enough, you can treat k<sub>1</sub> practically as μ. You can calculate for the new mean, lets call it μ[n+1] = (k<sub>n</sub>+Nμ[n]-μ[n])/N = μ[n+1] (k<sub>n</sub>)/N + ( (N-1)μ[n])/N If the signal ends up changing  from its original mean to a different value, the steady state of this method will converge to the newer mean. Our implementation is below:
-[Code Snippet]
+<code>
+function [Y] = RTprocess(X,N)
+%Given the input X which is rxk matrix the real time process output Y of
+%the peak detector of the past N samples and plots it.
+    %   where n is the moving average mean at sample n. It follows the
+    %   equation where nth element in vector Y is u[n], for each m columns
+    %   u[n] = (x[n]+u[n-1](N-1)u)/N. if x[n] doesn't vary too long and
+    %   N is large, it will take the approximate average of the last N
+    %   samples. If the signal is in steady state, this function will
+    %   converge to the average of the peak detector.
+X = abs(X-ones(size(X,1),1)*mean(X)); % Peak Detector x-> [-u]-> [|. |]-> y
+Y = zeros(size(X,1),size(X,2));       % Initialize output
+for n = 2:size(X,1)
+    Y(n,:) = X(n,:)/N + (N-1)*Y(n-1,:)/N; % u[n] = (x[n]+u[n-1](N-1)u)/N.
+    end
+plot(Y)                               % plot real time process
+title('Varied Input Trials')
+xlabel('sample')
+ylabel('system output')
+end</code>
+<code>
+figure(1);                                  % Plots simulational real time process
+RT = RTprocess(X,10000);
+legend([num2str(theta(1,:)') ones(size(theta,2),1)*(' degrees')])   %Attatch legend with units</code>
+{{:wind_sensor:rtangt.jpg|}}
 Notice if you don't chose a value N large enough, the output of the mean of the last N samples will very large, an example is shown below.
-[Code Snippet]
+<code>
-[Graph of non-tracked mean]
+% Test on how the N parameter makes a difference in find the moving average
-[Code Snippet]
+clc;clear;close all;                    % clears all
-[Graph of tracked mean]
+t = 0:1/26000:10;                       % Time vector
+x = 5*sqrt(t).*sin(10*t);               % signal with increasing mean
+figure(1)
+plot(t,x);                              % plot signal
+Label('Increasing Amplitude Signal','sec','Amplitude');
+figure(2)                               % Plot real time process with varied Ns
+RT = RTprocess(x',50000);
+title('real time process N = 50000')
+figure(3)
+RT = RTprocess(x',5000);
+title('real time process N = 5000')
+figure()
+RT = RTprocess(x',500);
+title('real time process N = 500')
+end</code>
+{{:wind_sensor:input.jpg|}}
+{{:wind_sensor:rt50000.jpg|}}
+{{:wind_sensor:rt5000.jpg|}}
+{{:wind_sensor:rt500.jpg|}}
 From this, we used the equation we found during the data correlation and fit it with the current output and estimate our wind speed. (Note: the model for features x<sub>1</sub>,x<sub>2</sub>,...x<sub>n</sub>, will output a model vector or [w<sub>1</sub> ,w<sub>2</sub> ,... w<sub>n</sub>] where y = (w<sub>1</sub>)(x<sub>1</sub>)+(w<sub>2</sub>)(x<sub>2</sub>)+...+(w<sub>n</sub>)(x<sub>n</sub>). If you use polyfit, with fitted N powers, the model would look like [w<sub>1</sub>, w<sub>2</sub>,...., w<sub>N</sub>, w<sub>N+1</sub>] and will be fitted as y = (w<sub>1</sub>)(x<sup>N</sup>)+(w<sub>2</sub>)(x<sup>N-1</sup>)+...+(w<sub>N</sub>)(x<sup>1</sup>)+(w<sub>N+1</sub>). A code below where we implement our code and how often we output is shown below.
@@ Line 75: / Line 161: @@
 Note: We did a piece wise regression at output<13 and output>13 because if you look above, that is when it changes a lot faster, and we rather avoid an extra power as a feature to decrease residual error and just use 2 models for 2 parts of the fitting.
-(Insert arduino code)
+The code for collecting data on Arduino is below. This is for collecting analog data before it is processed. This code collects data from all four microphones with a 1 microsecond delay in between readings. Note that each microphone is only being sampled one-fourth of the time. After releasing the reset button on the Arduino, microphone 1 will always be the first data point.
+<code>
+void setup() {
+  Serial.begin(14400);
+  analogReference(EXTERNAL);
+}
+void loop() {
+  delayMicroseconds(1);
+  int sensorValue1 = analogRead(A0);
+  delayMicroseconds(1);
+  int sensorValue2 = analogRead(A1);
+  delayMicroseconds(1);
+  int sensorValue3 = analogRead(A2);
+  delayMicroseconds(1);
+  int sensorValue4 = analogRead(A3);
+  Serial.println(sensorValue1);
+  Serial.println(sensorValue2);
+  Serial.println(sensorValue3);
+  Serial.println(sensorValue4);
+}
+end</code>