The purpose of this blog is to learn how to use the sasviya.ml and sklearn Python packages in SAS Viya Workbench to build and evaluate machine learning models. SAS Viya Workbench is a new SAS programming environment that supports the use of both the SAS and Python languages. The Workbench environment is designed to support native Python programming, and SAS has released a proprietary Python package named sasviya.ml that contains optimized SAS machine learning algorithms designed to run in SAS Viya Workbench. Of course, SAS Viya Workbench can also be used to execute SAS procedures, but this blog will focus on the Python functionality.   We created this example using a Python Jupyter notebook and the Visual Studio Code IDE, and the code in this example requires the use of SAS Viya Workbench. All of the code presented in this blog can be copied into a .py Python program or a .ipynb Jupyter Notebook in SAS Viya Workbench and you can execute it to follow along with the examples.   import requests # File path and name file_path = r"/workspaces/myfolder/MachineLearning/hmeq.csv" # Specify the URL of the CSV file url = r"https://support.sas.com/documentation/onlinedoc/viya/exampledatasets/hmeq.csv" # Download the and save CSV file to Workbench response = requests.get(url) with open(file_path, 'wb') as f: f.write(response.content) print(f'File downloaded:{file_path}')   To make this example more portable, we start by downloading data from a web URL using the Python requests package. We define a file location in the SAS Viya Workbench where we want to save the CSV, in this case in the /workspaces/myfolder/MachineLearning/ folder, and we want to save the file as hmeq.csv. We download the data from a SAS documentation page and use the requests package in Python to write the URL content to a CSV file. The /workspaces/myfolder/ location was created when we initialized our Workbench environment, and we just manually added the MachineLearning folder using the VSCode interface.   import pandas as pd hmeq_df = pd.read_csv(r"/workspaces/myfolder/MachineLearning/hmeq.csv")   Although we will be using SAS Viya Workbench machine learning algorithms, this example is all in Python, so we will use Pandas and Scikit-Learn to help load and prepare the data for machine learning. The SAS machine learning algorithms all live in the sasviya.ml Python package, and we will import them individually as we use them. We import pandas and read the CSV file into a pandas dataFrame named hmeq_df. SAS Viya Workbench includes commonly used Python packages for data science, and in general you can use PIP to install your own packages, although the write permissions and available packages will be determined by your site administrator. In this example we use pandas, numpy, and scikit-learn to prepare and evaluate our machine learning models. We didn’t have to pip install these packages; they were already included as part of the configuration of SAS Viya Workbench.   Select any image to see a larger version. Mobile users: To view the images, select the "Full" version at the bottom of the page.   Before we start building models, let’s quickly build a little bit of intuition about this data. This is fictitious financial services data; we operate a bank and must decide whether to accept or reject loan requests from customers. Our target in this data is the variable BAD, or binary applicant default and it indicates whether our historical customers defaulted on their loan or not. Here BAD equals 1 indicates that the customer did default on the loan; BAD equals 0 indicates they did not default on the loan. The rest of the variables in the data represent information about the customer collected by the bank, which can be used as candidate input variables for our machine learning models. These input columns contain information like the amount of the home equity loan the customer is requesting, the amount left on their mortgage, the value of the mortgaged home, the reason they are requesting a loan, their job, their years on the job, along with some information about their past credit history like the numbers of derogatory credit reports and delinquent credit lines. The important thing about this data from a technical software perspective is that we have a mix of interval variables like LOAN and MORTDUE along with categorical variables like REASON and JOB (in our data these variables contain strings). We don’t have to do any dummy coding when using the sasviya.ml package, but the SAS machine learning models will treat the categorical input variables different from interval variables.   hmeq_df.head(10)     We use the head() method to print out the first 10 observations in our dataset. Looking at these 10 printed rows, we can see that we have some missing values in the data, in Python coded as numpy NaN values. We want to do something to replace these missing values, otherwise the rows with missing values won’t be included in our machine learning models, which can dramatically reduce the amount of data we have available for training the models. Imputing these missing values will be a part of our data preprocessing, which is our next step. This is the most basic preprocessing that we must do for machine learning (the bare minimum). There are lots of data preprocessing methods available in Python packages in SAS Viya Workbench that could be useful with different kinds of data, and in general when you are working with new data it is recommended to spend some time on data exploration, which can help better identify what kind of data preprocessing would be useful. For this example we will stick with simple data preprocessing, which will include partitioning the data, imputing missing values, and selecting relevant inputs variables. Partitioning the data into training and validation samples is a best practice that you cannot skip in machine learning; if you don’t have a validation sample you can’t really know if you have fit an effective model or if the model has just memorized the historical data used for training. We impute missing values; we do have a lot of missing values in the data, and we want to use as much of the data as possible. Automatic variable selection is more useful when you have a lot of hard-to-understand input variables (we only have 12 easy-to-understand input variables in this dataset), but we include it in this example because it can be useful to compare your intuition about what variables might be useful in predicting the target (manual variable selection) with the results of data preprocessing routines that automatically select variables.   #We use sklearn to partition the data into training and validation data from sklearn.model_selection import train_test_split X = hmeq_df.drop('BAD', axis=1) y = hmeq_df['BAD'] X_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.3, random_state=919) print("Training data shape:", X_train.shape, "\n" + "Validation data shape:", X_valid.shape)   We split the HMEQ dataset into training and validation samples using the train_test_split function from the scikit-learn model selection package, sampling 30% of the data as a validation set and the remaining 70% as a training set. Scikit-learn expects us to split the inputs from the target, so we create a dataFrame X containing only the inputs and a dataFrame y containing the target BAD. We apply the train_test_split function to these dataFrames, specifying that we want a 30% test sample for the validation data, and that we want to stratify based on the target variable. This stratified sampling means that we will have the same proportion of BAD=1 (people who defaulted on their loans) in the training and validation data. These samples are supposed to represent the data our machine learning model will encounter in deployment, so we should preserve the target distribution from the original dataset. We also specify a random state to ensure that the sampling is reproducible when we run this notebook multiple times. This process outputs 4 different pandas objects, 2 dataFrames containing training and validation inputs (x_train and x_test), and 2 series containing training and validation targets (y_train and y_test). After the partitioning we have 4,172 training observations and 1,788 validation observations.   Now that we have built this partition, we can start using the training data to build models and do our data preprocessing, while remembering to apply all the data preprocessing steps to the validation data. Our first form of data preprocessing is to impute missing values, but we want to use the training data to calculate the median and the mode that we're going to use for imputation, and then we're going to use those values when imputing the validation data. The idea is that it's cheating to use information from the validation data to impute values on the training data, since information leaks from the validation sample into the training process. We need to learn our parameters from the training data, but the median and mode that we use for imputation can be thought of as model parameters associated with data preprocessing. This is also more reflective of a deployment scenario, where we will also use the median and mode from the training data to impute missing values in new observations that we want to score.   #We use sklearn to impute missing values (using the median/mode from the training data to impute the validation data) from sklearn.impute import SimpleImputer import numpy as np #Calculate the median and mode on the training data to prepare the imputer imp_interval = SimpleImputer(missing_values=np.nan, strategy='median') imp_interval.set_output(transform='pandas') imp_interval = imp_interval.fit(X_train.select_dtypes(include='number')) imp_nominal = SimpleImputer(missing_values=np.nan, strategy='most_frequent') imp_nominal.set_output(transform='pandas') imp_nominal = imp_nominal.fit(X_train.select_dtypes(include='object')) #Apply the imputation to the training data X_train_imp_int = imp_interval.transform(X_train.select_dtypes(include='number')) X_train_imp_nom = imp_nominal.transform(X_train.select_dtypes(include='object')) X_train_imp = pd.concat([X_train_imp_int, X_train_imp_nom], axis=1) #Apply the imputation to the validation data X_valid_imp_int = imp_interval.transform(X_valid.select_dtypes(include='number')) X_valid_imp_nom = imp_nominal.transform(X_valid.select_dtypes(include='object')) X_valid_imp = pd.concat([X_valid_imp_int, X_valid_imp_nom], axis=1) print("Training data shape:", X_train_imp.shape, "\n" + "Validation data shape:", X_valid_imp.shape)   We import the SimpleImputer function from the sklearn.impute package, and we create two imputer objects, one for the interval variables that imputes using the median, and one for the nominal variables that imputes using the mode (‘most frequent’ is the term used by sklearn). For both imputers we set the output to be a Pandas dataFrame (which is not the default for the imputers), this will be important later because the sasviya.ml algorithms are easier to use when our data is in the form of a Pandas dataFrame. In both cases we use the fit method on the training data to learn the imputation parameters (in this case just the median and mode for the columns with missing values) based only on the training data. It is important to separate the interval and the nominal variables for this, so we use the select_dtypes method to ensure that the median imputation is calculated for the numeric variables and the mode imputation is calculated for the textual variables. This could be a bit more complicated if you have nominal variables with numeric values, at which point it could be helpful to create lists of input variables separated by how they should be treated. Fitting the two imputer objects that we created just learns the model parameters from the training data, but we also must apply the imputation to both the training and validation data. We do this separately for the training and validation data, for each sample concatenating the interval and nominal input variables into a single dataFrame. At the end of this process, we still have the same size training and validation samples, but this time the 12 input columns have no missing values. We also have the imputer objects that we can use to apply the imputation to new data during the scoring process.   #Select useful input variables using chi-square test from sklearn.feature_selection import SelectKBest, chi2 #Sklearn requires nominal variables to be dummy coded, for now let's skips this and just select from the interval variables #Select variables based on the training data selector = SelectKBest(chi2, k=8).fit(X_train_imp_int, y_train) selector.set_output(transform="pandas") #Apply the selection to the training and validation data train_selected = selector.transform(X_train_imp_int) valid_selected = selector.transform(X_valid_imp_int) X_train_final = pd.concat([train_selected, X_train_imp_nom], axis=1) X_valid_final = pd.concat([valid_selected, X_valid_imp_nom], axis=1) print("Training data shape:", X_train_final.shape, "\n" + "Validation data shape:", X_valid_final.shape) print("Selected interval features:", selector.get_feature_names_out()) print("Selected nominal features:", list(X_train_imp_nom.columns))     Variable selection is important for ensuring that the machine learning models are provided useful inputs. You can always select important variables based on subject matter expertise, but it’s often helpful to compare your intuition to an automatic variable selection approach. We import the SelectKBest function from scikit-learn to select the top 8 most relevant interval variables, based on a chi-square test. We only have 10 interval input variables so this will drop the 2 interval variables that are least useful in predicting the target. Scikit-learn requires categorical variables to be dummy coded (often using one-hot encoding) before they can be used for most machine learning and data preprocessing, but the sasviya.ml package allows us to skip this step so we will leave our categorical variables in their original form. For this example, we will just keep the 2 categorical variables as inputs and select the 8 best interval inputs using SelectKBest. Just like before we want to output a pandas dataFrame, and once we select variables based on the training data (once we fit the SelectKBest method), we apply the variable selection to the interval variables in both the training and validation data. After reducing the interval variables, we concatenate them back with the categorical variables to get the preprocessed training and validation input samples. This time these samples only have 10 columns, since we have eliminated 2 of the variables. We also print out the selected variables, but we didn’t apply any variable selection to the categorical variables, so we selected both without doing any calculations. We can now use these preprocessed input samples to fit our machine learning models (with the training data) and evaluate and compare the models’ performance (with the validation data).   The sasviya.ml package is designed to integrate with scikit-learn, so the sasviya.ml model objects have the same syntax and most of the same functionality as the equivalent scikit-learn model objects. The major difference is in the execution, the models in the sasviya.ml package execute using optimized SAS libraries, which takes advantage of multithreading and can yield much faster runtimes than the equivalent models in scikit-learn. If you are interested in comparing the performance of sklearn and sasviya.ml, there is a link to a repository of speed comparisons in the references at the end of this blog. The sasviya.ml syntax is nearly identical to the scikit-learn syntax (the constructor is the same, the objects have the same methods such as fit and transform, etc.), so if you have experience building machine learning models in scikit-learn, you can easily start building models using sasviya.ml. One convenient feature about this design is that you can take existing scikit-learn code and without modifying the actual code, you can change the import statement to point to sasviya.ml instead of scikit-learn. This will use the optimized SAS libraries for machine learning without requiring you to do any code rewrites. One difference between scikit-learn and the sasviya.ml package to keep in mind is that the scikit-learn models require all inputs to be numeric (so we must dummy code before sending the data to the model), whereas the sasviya.ml models are designed to do the dummy coding for you if you provide categorical inputs. This is a convenience feature, so all your legacy scikit-learn code will work with sasviya.ml, but you can skip a few steps in code when using the SAS optimized libraries. The rest of this notebook will look just like a scikit-learn demo, but we will import the models from sasviya.ml instead of from sklearn. If you have experience with scikit-learn and would prefer to organize your code differently or use features like sklearn pipelines, you can easily import your sklearn code into SAS Viya Workbench and modify it to use the optimized SAS algorithms.   #fit a simple logistic regression model from sasviya.ml.linear_model import LogisticRegression logreg = LogisticRegression(solver='lbfgs', tol=1e-4, max_iter=1000) logreg.fit(X_train_final, y_train)   Let’s start building some different machine learning models, always making sure to compare to a simple logistic regression model. Logistic regression isn’t exactly what we think of when we think of a machine learning model, but it’s always a good idea to compare our complex nonlinear models to a simple linear model. We import the LogisticRegression object from sasviya.ml.linear_model and use the constructor to instantiate a logistic regression object named ‘logreg’. We accept most of the default settings for the logistic regression model, specifying that we want to use the LBFGS method for the likelihood maximization, with a convergence tolerance of 10 -4 and a maximum of 1000 iterations for the LBFGS solver. These settings are all about how the model estimates the parameters, but there is an option for choosing a variable selection method (like backward or stepwise) for the logistic regression model. After we construct the ‘logreg’ object, we can apply the fit method to the training data (X_train_final and y_train) to create our model. After the fitting process the ‘logreg’ object can be used with the transform method to score new data (which we will do later in the notebook when we evaluate model performance on the validation data).   #fit a decision tree model from sasviya.ml.tree import DecisionTreeClassifier dtree = DecisionTreeClassifier(criterion="chisquare", max_depth=10, ccp_alpha=0) dtree.fit(X_train_final, y_train)   We can do the same thing with a decision tree model, importing the DecisionTreeClassifier object from sasviya.ml.tree. We use the constructor to instantiate the model object ‘dtree’, and this time we specify some key decision tree hyperparameters. We use the chi-square method as the split criterion, which is how we evaluate the quality of the splits in the decision tree. We also specify a max depth of 10, so any leaf node in the tree is no more than 10 splits below the root node. A smaller max depth leads to a simpler tree, which can be useful to avoid overfitting, but can also lead to a model that is too simple to represent the training data. Just like with the logistic regression model, we fit the dtree object using the training data, and we will use this same object later to evaluate performance on validation data.   #fit a random forest model from sasviya.ml.tree import ForestClassifier forest = ForestClassifier(criterion="chisquare", n_estimators=100, max_depth=7, min_samples_leaf=5, bootstrap=0.6, random_state=919) forest.fit(X_train_final, y_train)   Next, we create a random forest model by importing the ForestClassifier object from sasviya.ml.tree. A forest model is an ensemble of individual decision trees, with each tree fit to be slightly different from the other trees in the forest. When we create the forest model, we specify some hyperparameters related to how individual trees are grown, like max_depth, criterion, and min_samples_leaf. These hyperparameters are the same as the ones we use with the individual decision tree model.  We can also specify hyperparameters associated with how we create the forest, like n_estimators and bootstrap. The forest model is only useful if the trees in the ensemble are different from one another, so we set the bootstrap options to 0.6 so that each tree is fit using a random sample with 60% of the training data. Once we use the constructor to define the model hyperparameters, we can fit the forest object using the training data.   #fit a tree-based gradient boosting model from sasviya.ml.tree import GradientBoostingClassifier gradboost = GradientBoostingClassifier(n_estimators=100, max_depth=4, min_samples_leaf=5, learning_rate=0.1, subsample=0.8, random_state=919) gradboost.fit(X_train_final, y_train)   We can also fit a gradient boosting ensemble model, where instead of fitting a forest of independent trees, we fit a sequence of trees where each tree tries to improve on the performance of the previous trees. The constructor for the GradientBoostingClassifier object has many of the same hyperparameters as the forest and decision tree models, along with some options specific to creating the sequence of trees like the learning_rate hyperparameter. Notice with this model we use a much smaller maximum depth than with the previous models; in general, gradient boosting models are often more powerful and accurate than individual trees, but this comes with the risk of overfitting the training data. We use a small max depth for this model to help ensure that we don’t memorize the training data. Once we construct the gradboost object, we fit it using the training data.   #fit a support vector machine classifier from sasviya.ml.svm import SVC svm = SVC(C=1.0, kernel="rbf") svm.fit(X_train_final, y_train)   Finally, we fit a support vector machine model (in this case a support vector classifier, or SVC model which we import from sasviya.ml.svm). For binary targets this is normally a linear model that finds a hyperplane to maximize the separation between the two classes. In this example we use  the radial basis kernel function (known as the RBF kernel), which is a nonlinear kernel that projects the problem into a higher-dimensional space and then fits a linear hyperplane in this projection. In the original dimensions this becomes a nonlinear decision boundary. In addition to the kernel, we specify the penalty parameter C, which determines a cost penalty for incorrectly classified cases. Once we create the svm model object, we fit it using the training data.   Now that we have fit our five different machine learning models using training data (well really four machine learning models and a logistic regression as a baseline comparison) we can evaluate the performance on validation data and choose a champion model. The example we show in this blog uses a binary target, but each of the models we used also has an option for data with a continuous target. For most of the models we have a regressor option, so instead of using the ForestClassifier model, we would use the ForestRegressor model for a continuous target. More details about these regressor models are available in the SAS Viya Workbench documentation.   #score training and validation data using the fitted models (we can assess them separately after this) models = ['logreg', 'dtree','forest','gradboost','svm'] train_out = dict.fromkeys(models, None) valid_out = dict.fromkeys(models, None) #score the training and validation data using the models, and join the predictions to the target for model in models: train_out[model] = eval(model).predict_proba(X_train_final) train_out[model] = train_out[model].join(y_train) valid_out[model] = eval(model).predict_proba(X_valid_final) valid_out[model] = valid_out[model].join(y_valid)   To evaluate model performance and ensure we don’t overfit the training data, we score both the training and validation data using all our models. First, we create a list of models, with the name of the model on the list equal to the name of the model object we created during training. We then create an output dictionary for training and validation data from these models, so we have a place to store the scored training and validation data. We use a Python for loop to loop over the list of models, generating the predicted probabilities using the predict_proba method applied to the model object. We join these predicted probabilities with the true value of the target and save the scored output in the model dictionary. We do this separately for training and validation data to yield two output dictionaries, each containing dataFrames with scored output for each model.   #plot ROC from sklearn.metrics import roc_curve from matplotlib import pyplot as plt train_roc = dict.fromkeys(models, None) valid_roc = dict.fromkeys(models, None) plt.figure(figsize=(10,8)) plt.xlabel('False Positive Rate') plt.ylabel('True Positive Rate') plt.title('ROC Curve for Training Data') for model in models: train_roc[model] = roc_curve(train_out[model]['BAD'], train_out[model]['P_BAD1'], pos_label=1) plt.plot(train_roc[model][0], train_roc[model][1]) plt.legend(models) plt.figure(figsize=(10,8)) plt.xlabel('False Positive Rate') plt.ylabel('True Positive Rate') plt.title('ROC Curve for Validation Data') for model in models: valid_roc[model] = roc_curve(valid_out[model]['BAD'], valid_out[model]['P_BAD1'], pos_label=1) plt.plot(valid_roc[model][0], valid_roc[model][1]) plt.legend(models);       We plot the ROC curve by using the roc_curve function in the sklearn.metrics package to calculate ROC information. Once again, we create empty dictionaries to store ROC information, and then we loop over each of our models, adding the ROC curve information for each model to the pyplot object. We do this separately for training and validation data so we can compare models on both partitions. There are plenty of other assessments we could use, and we could also use other Python plotting packages like seaborn. Looking at the results we can see that the gradient boosting model is the champion on both training and validation data, and although the performance is better on the training data, we don’t see evidence of major overfitting, the validation and training curves seem similar in the plots.   #Calculate misclassification for 'champion' model on validation data from sklearn.metrics import accuracy_score cutoff = 0.5 valid_preds = [1 if valid_out['gradboost']['P_BAD1'][elem] > cutoff else 0 for elem in valid_out['gradboost'].index] print('Misclassification Rate on Validation Data for Gradient Boosting at',cutoff,'cutoff:', 1-accuracy_score(valid_out['gradboost']['BAD'], valid_preds))       Before we declare the gradient boosting model the final champion model, let’s look at the misclassification rate at the 0.5 probability cutoff. This isn’t necessarily the best cutoff for all models on all data, but it is a good place to start when building and evaluating models. We import the accuracy_score function from sklearn.metrics, but before we can apply it to the output classification scores, we have to apply the 0.5 probability cutoff to the predicted probabilities generated earlier using the predict_proba() method. We could have used the predict() method to generate these class values instead, but it’s helpful in this example to explicitly choose a cutoff (mainly so we can emphasize that the cutoff affects the outcome and has an impact on the misclassification rate). We use Python list comprehension to create a list of class predictions on the validation data, and then we apply the accuracy_score function to this list, comparing it to true value of the target. Of course, we print the misclassification rather than the accuracy, so we see that the champion gradient boosting model is incorrect for about 9.34% of the observations in the validation data.   The last step is to start thinking about model deployment. To deploy these models, we just need to apply any data preprocessing and the champion model to the new data. It’s easy to score new cases in the same SAS Viya Workbench environment, we just keep the trained data-preprocessing/model objects in Python memory and apply the transform/predict methods to the new data. It gets a bit more complex if we want to deploy to a different environment; in this case we would want to pickle the Python model objects and then unpack them in a different environment.   References: SAS Viya Workbench Documentation SAS Viya Workbench Python Programming Guide ... View more

The purpose of this blog is to learn how wavelets can be used to remove unwanted high-frequency noise from digital signals. Wavelets can be used to decompose a signal into different frequency components (with different resolutions at each frequency level), and once we have the decomposition, we can manipulate the high-frequency component of the signal to reduce or eliminate the noise. This can be useful to clean up sensor data with unexpected high frequency noise. An example of this is a vibration sensor on a spinning turbine that generates an electromagnetic signal, but the signal experiences electromagnetic interference before it is read into a computer. The output digital signal includes both the actual vibration data, and some unwanted electromagnetic noise riding on top of the signal.   The basic process for this is to calculate the wavelet decomposition coefficients, and then apply some kind of threshold to the wavelet detail coefficients (the coefficients for the high-frequency components of the signal) to shrink their magnitude (soft thresholding) or reduce them to zero (hard thresholding). We will describe this process in detail in this blog and we will see examples of applying this to simulated and real data using the wavelet tools in the SAS Interactive Matrix Language.   In a previous blog on wavelet analysis (link provided below), we introduced the discrete wavelet transform and showed the following equation for a family of wavelets at different shift and scale parameters (scaling by powers of two):     The scale parameter j determines the frequency of the wavelet, while the shift parameter k determines the temporal support for the wavelet. As we increase k, we are moving the wavelet forward in time (to the right on a plot) across the input signal. As we increase j, the frequency of the wavelet increases and thus we capture higher frequency variation in the signal using wavelets with larger values of j. The wavelet transform yields a collection of coefficients representing the information in the signal corresponding at different shifts and scales, with one coefficient for each wavelet (each pair of values of k and j define a wavelet and thus correspond to one coefficient). The coefficients are thus spread out in time and frequency, with different values of k corresponding to the signal components at different times, and different values of j corresponding to the signal frequency contents. The coefficients for high values of j are called the “detail coefficients” and represent the high-frequency components in the signal. We will threshold these detail coefficients to eliminate some of the high-frequency components and thus remove unwanted high-frequency noise.   Let’s explore an example of removing low-amplitude, high-frequency noise from a simulated sinusoidal signal. This is a very simple toy example (most real-world signals are not perfect sinusoids with perfect sinusoidal noise) but it provides good intuition for how the denoising operation works. We start by simulating the input signal and the additive noise using SAS/IML:   proc iml; /*We start by simulating a digital signal with a small amplitude high-frequency noise injected*/ f_s = 500; /*sampling frequency*/ f_s2 = f_s/2; /*Nyquist frequency*/ t_s = 1/f_s; /*sampling period (time between samples)*/ N_samp = 1000; /*number of samples*/ t = do(0, (N_samp-1)*t_s, t_s); fS = 2; fN = 225; pi = constant("pi"); x_Signal = sin(2*pi*fS*t); /*Signal component*/ x_Noise = 0.05*sin(2*pi*fN*t); /*Noise component*/ x_NoisySignal = x_Signal+x_Noise; create work.input_signal var {"t","x_NoisySignal","x_Signal","x_Noise"}; append; close work.input_signal; /*plot our signal in the time domain*/ submit; title "Input Digital Signal with and without Noise"; proc sgplot data=work.input_signal; series x=t y=x_NoisySignal; series x=t y=x_Signal; xaxis grid label="Time (s)"; yaxis grid label="Amplitude"; legend; run; endsubmit; run;   Note that in this example the input signal has a frequency of 2 Hz, while the additive noise has a frequency of 225 Hz. The noise is also much smaller than the signal in this example, with an amplitude of 5% the signal amplitude. This yields the following plot, showing the noise-free input signal overlaid on top of the actual noisy signal we will analyze:       Our goal is to reconstruct the original noise-free signal (x_Signal in the plot) from the ‘measured’ input signal with noise (x_NoisySignal in the plot). In this example we are imagining we don’t have access to the noise-free signal, and the only thing we have is the measured data, x_NoisySignal. We start by performing the discrete wavelet transform to calculate the wavelet decomposition, and then we plot the detail coefficients so we can explore them before we apply any kind of thresholding (note that this code is still in the same call to PROC IML as the previous block of code. Listed macro and function calls are described below):   /*now let's calculate the wavelet decomposition so we can threshold the wavelet detail coefficients*/ %wavinit; /*let's use a predefined Daubechies 3-tap wavelet*/ optn = &daubechies3; /*now perform the discrete wavelet transform to find the decomposition*/ call wavft(decomp, x_NoisySignal, optn); call wavprint(decomp, &detailCoeffs, 1, 4); call coefficientPlot(decomp); call coefficientPlot(decomp, &SureShrink); run;   The %wavinit macro is used to setup the wavelet functions we will use to streamline the submission of the wavelet analysis code. This code allows us to define our wavelet family and type using the macro variable &daubechies3 (as opposed the syntax we would use without the %wavinit macro which would require looking up all of the numeric options in the SAS documentation and specifying each option by number). The function call to wavft calculates the wavelet decomposition, but we will use other wavelet analysis functions to extract the relevant information from this decomposition. We start by printing the detail coefficients at level 1-4 using the wavprint call function. This yields the following table of coefficients:       The different value of the Translate column correspond to different shift values, but they fill the time domain for each level. Although the values of Translate represent the shift k, they don’t correspond to the actual shift value in time. In this way the two detail coefficients at level 1 cover the same time domain as the four detail coefficients at level 2, we just have greater time resolution at the higher levels. This table only shows the first four levels of detail coefficients, but with our default settings we calculate 10 levels of detail coefficients (levels 0 to 9). We want to apply a threshold to the higher-level detail coefficients, but there are too many of them to explore in a table, so we plot them:       This plot clearly illustrates how the detail coefficients are distributed in time (with greater time resolution for the more numerous high frequency detail coefficients). We haven’t applied thresholding yet, so we can see the high-frequency variation in the higher levels of detail coefficients. This plot also highlights the time-frequency resolution discussed in the previous blog on wavelet scalograms, with greater time resolution at higher frequencies (although it doesn’t highlight the corollary to this, which is that we have lower frequency resolution at higher frequencies as well).   We can apply either soft thresholding to these coefficients to shrink them in magnitude, or hard thresholding to set some of the coefficients to exactly zero. We use the &SureShrink option in SAS/IML to apply soft thresholding to these detail coefficients, although we could also try out the &RiskShrink option to apply hard thresholding. The software allows us to use more customized thresholding options, but the developers include macro variables in the wavelet toolkit to make it easier to apply popular thresholding methods. After applying the soft thresholding, we can plot the transformed detail coefficients:     It looks like the detail coefficients at levels 6, 7, 8, and 9 have been set to zero, although we used a soft threshold, so they have actually just been shrunk to the point where they are not visible on the plot. We can see that most of the high-frequency variation we saw in the detail coefficients at the higher levels is absent from the plot, so we can expect that this thresholding has removed most of the high-frequency noise. To investigate this, we must perform the inverse wavelet transform on the decomposition with the modified detail coefficients:   call wavift(x_denoised, decomp, &SureShrink); x_diff = x_signal - x_denoised; /*create an output dataset to plot the results*/ create work.denoised_signal var {"t","x_NoisySignal","x_Noise","x_Signal", "x_denoised", "x_diff"}; append; close work.denoised_signal; submit; /*first we plot the denoised signal*/ title "Denoised Signal"; proc sgplot data=work.denoised_signal; series x=t y=x_denoised; xaxis grid label="Time (s)"; yaxis grid label="Amplitude"; run; /*now we plot the difference between the original "noise-free" signal and the denoised signal*/ title "Difference between original noise-free signal and denoised signal"; proc sgplot data=work.denoised_signal; series x=t y=x_diff; xaxis grid label="Time (s)"; yaxis min=-1 max=1 grid label="Amplitude"; run; endsubmit; quit;   The wavift call performs the inverse wavelet transform to create the output signal x_denoised from the wavelet decomposition decomp, using the &SureShrink method to threshold the wavelet detail coefficients. After we create the denoised signal, we calculate a difference between the original noise free signal and the denoised signal to evaluate the effectiveness of the wavelet denoising approach. We plot the denoised signal along with the difference between the denoised signal and the original noise-free signal to compare:     Visually the denoised signal looks like our original signal; it no longer has the high-frequency variation we saw on the noisy signal. Of course, the plot can be misleading, so we can also plot the difference between the denoised signal and the original noise-free signal:     We can see that the denoised signal is not exactly equal to the original noise-free signal, and the deviations seem to be somewhat periodic in time (which makes sense given the periodic nature of the wavelet analyses).   This toy example is just a simple sinusoid, but wavelet denoising is especially useful when the data contains sharp transitions or jumps, something that is very common in real world sensor and signal data. Traditional Fourier transform methods can also be used to remove high-frequency noise (discussed in a previous blog post, link below), but the basis used for the Fourier transform (sines and cosines) does a bad job representing sharp transitions or level shifts. Any attempt to use a Fourier transform to remove high-frequency noise will also invariably smooth out and modify these sharp transitions. This is a problem because these sharp transitions are not high-frequency noise components, and they often represent times when important events occurred in the data, so we want to preserve the temporal location of these transitions. Wavelet analysis is a powerful tool in these scenarios because the wavelet transform can remove the high-frequency noise without smoothing or eliminating regions where the signal experiences a level shift or a sharp transition.   References: Previous blog post on time-frequency analysis using Short-Time Fourier Transform Previous blog post on time-frequency analysis using Wavelet Transform SAS Documentation describing wavelet analysis SAS Documentation page listing the wavelet options and symbolic names ... View more

The purpose of this blog is to learn how wavelets can be used to analyze signal and sensor data. We learn what wavelets are and how they can be used to decompose periodic signals into their constituent frequencies. In this blog we focus on the basics of wavelet analysis and compare the time-frequency decomposition generated by wavelet analysis (we will explore a plot of this decomposition called the ‘scalogram’) with the spectrogram generated by traditional Fourier analysis (using the short-time Fourier transform). This is just one application of wavelets, with other popular applications including denoising and data compression. We will explore other wavelet applications in subsequent blog posts.   Wavelets are short oscillations with a finite duration in time. There are a variety of families of functions that serve as wavelets, with the common theme being a function that starts at zero and ends at zero with some oscillations in between. This is in contrast with the basis functions used in Fourier analysis, the sine and cosine waves that have infinite extent. Thus, from the perspective of Fourier analysis, wavelets are already “windowed” in that they are zero valued outside of a small, finite duration. Plotted below are two examples of simple wavelets (more complex wavelets are often used in practice). The first plotted wavelet is a Morlet wavelet and the second is a Haar wavelet.    Morlet Wavelet Select any image to see a larger version. Mobile users: To view the images, select the "Full" version at the bottom of the page.      Haar Wavelet   These “mother wavelets” are used to define a collection of similar wavelets at different scales and locations in time. From a given mother wavelet, we derive wavelets at different shifts and scales. Correlations between these wavelets and our signal of interest identify the signal components at different frequencies (scales) and at different points in time (shifts). This idea of ‘scaling’ and ‘shifting’ the wavelet is key to calculating the continuous wavelet transform (although for digital signals we will have to calculate the discrete wavelet transform), which is the time-frequency decomposition of our signal. Given a mother wavelet, ψ, we can define the child wavelet at scale a and shift b as follows:     We can think of the scale a as representing frequency; as a increases, the frequency of the oscillation in the wavelet increases. We can think of the shift b as representing time; as b changes, we move the wavelet forwards or backwards in time relative to the original signal. Of course, this all applies to continuous wavelets. In the discrete case (anytime we are dealing with digital signals we must use the discrete wavelet transform), we define the child wavelets based on integers m and n as follows:     In this case the integer m is associated with the scale and the integer n is associated with the shift. The most common implementation of the discrete wavelet transform involves shifting and scaling by powers of two, enabling much more efficient algorithms for computing the transform. We thus can simplify the equation above:     In this case the integer j is the scale parameter, and the integer k is the shift parameter. The set of child wavelets for integer values of j and k establish a complete orthonormal basis for square-integrable functions of a real variable. In other words, they form a basis for the kinds of functions we encounter when working with signal data. This is analogous to how sines and cosines form a complete basis when working with the Fourier transform. The basis used in wavelet analysis can do a better job balancing time and frequency resolution when decomposing our signal as compared to the Fourier transform basis. We’ll talk about why this is in more detail below.   When we perform the short-time Fourier transform, window functions are applied to break the input signal into a bunch of different, finite duration portions. The Fourier transform is then applied to each portion to calculate the frequency spectrum. Since each portion represents a different time point in the original signal, we can assemble the frequency spectra into a two-dimensional spectrogram where we plot time on the x-axis (corresponding to the different windowed portions of the signal) and we plot frequency on the y-axis (the information that was contained in each frequency spectrum). The time resolution we plot on the x-axis are determined by the location and number of time windows we apply to the input signal while the frequency resolution is determined by the length of signal in the windows (the size of the windows). Although we can seemingly improve frequency resolution in the Fourier transform by zero-padding, this just interpolates values in the frequency spectrum; the fundamental frequency resolution of the transform is determined by the length of the signal. Thus, there is a fundamental tradeoff in time-frequency resolution in the short-time Fourier transform. Longer windows mean more frequency resolution, but with a fixed-length signal longer windows means we have fewer windows overall, and thus the time resolution is limited. Shorter windows allow us to improve the time resolution by including more windows to cover the signal, but of course degrade the frequency resolution. This fundamental limitation is the same basic idea as the Heisenberg uncertainty principle in quantum mechanics (we cannot simultaneously identify the exact location and momentum of a particle) and is often called the Gabor limit in signal processing. Wavelet analysis avoids some of the consequences of this limitation by enabling us to vary the time-frequency resolution at different locations in the spectrogram (which we will now call a scalogram).   To calculate the wavelet scalogram for a digital signal, we find wavelet coefficients for the signal associated with each of the shifted and scaled wavelets we derived previously from the mother wavelet. These coefficients are analogous to the power in the time-frequency bins from the STFT and can be calculated by finding the dot-product between the input signal and the corresponding child wavelet (at a specific value of the shift and scale parameters). In this way we don’t have to window the input signal at all, since each of the child wavelets are already finite in duration (as opposed to the sines and cosines used in the STFT). The magnitude of these wavelet coefficients are plotted in a two-dimensional scalogram with the shift parameter representing the different times and the scale parameter representing different frequencies. The key advantage over the STFT results from the construction of the child wavelets from the mother wavelet. We have longer wavelets at lower frequencies (thus we have bad time resolution but good frequency resolution) and we have shorter wavelets at higher frequencies (thus we have better time resolution but worse frequency resolution). This allows us to more easily tell when unwanted high-frequency noise is included in a signal, while still allowing us to accurately identify long-lasting low-frequency oscillations in the signal. A toy visualization of the time-frequency resolution tradeoff in the STFT versus the wavelet transform is displayed in the following figures:     In the STFT the time-frequency resolution is fixed across all frequency scales.     In the wavelet transform we can have longer wavelets at lower frequencies and shorter wavelets at higher frequencies, allowing the time-frequency resolution to vary across frequency scales.   Let’s explore some examples of plotting time-frequency decompositions in SAS/IML. We will start by simulating a digital signal with a linearly increasing frequency and then plot the Fourier transform spectrogram and the wavelet transform scalogram.   proc iml; /*start by simulating the chirped signal*/ f_s = 1000; /*sampling frequency*/ f_s2 = f_s/2; /*Nyquist frequency*/ t_s = 1/f_s; /*time between samples*/ N_samp = 2000; /*number of samples*/ pi = constant("pi"); t = do(0, (N_samp-1)*t_s, t_s); /*array of times*/ /*create the chirp*/ f0 = 10; /*starting frequency*/ c = 50; /*chirp, i.e. rate of change of the frequency*/ f_t = f0 + c*t; /*here is the linear chirp*/ /*f_t is the instantaneous frequency*/ phi_t = 2*pi*(f0*t + c/2*t#t); /*we integrate the linear chirp w.r.t. time and multiply by 2*pi to get the phase*/ /*phi_t is the instantaneous phase*/ x_in = sin(phi_t); create work.input_signal var {"t","x_in","f_t"}; append; close work.input_signal; /*plot the instanaeous frequency and our signal in the time domain*/ submit; title "Instantaneous Signal Frequency"; proc sgplot data=work.input_signal; series x=t y=f_t; xaxis grid label="Time (s)"; yaxis grid label="Frequency"; run; title "Input Signal with Linear Chirp"; proc sgplot data=work.input_signal; series x=t y=x_in; xaxis grid label="Time (s)"; yaxis grid label="Amplitude"; run; endsubmit; /*now we do a short time fourier transform on the data to explore how the frequency changes*/ winLen = 500; /*this will have a big impact on the time/frequency resolution*/ window = tfwindow(winLen, "Gaussian"); noverlap = floor(0.6*winLen); nfft = 512; /*increasing this increases the number of frequency bins*/ ramp = palette("Spectral",7)[7:1]; /* reverse "spectral" color ramp */ title "Short-Time Fourier Transform Spectrogram"; call spectrogram(x_in,window,noverlap,nfft) scale="linear" colorramp=ramp sampRate=f_s; /*next we do a wavelet transform and plot a scalogram /*the wavelet module expects a column vector*/ x_in = x_in`; /*initialize the wavelet graphics module*/ %wavinit; /*set wavelet options using macro variables from %wavinit macro*/ /*in this case let's use a predefined Daubechies 3-tap wavelet*/ optn = &daubechies3; /*now perform the discrete wavelet transform to find the decomposition*/ call wavft(decomp, x_in, optn); /*plot a scalogram of the wavelet decomposition*/ title "Wavelet Transform Scalogram"; call scalogram(decomp); quit;   This yields the following spectrogram and scalogram plots:       In this example the spectrogram does a better job of revealing the simple linear trend in the frequency, but we can still see the linearly increasing pattern the wavelet coefficients. Notice how the highest frequency detail coefficients in the wavelet scalogram have a high time resolution but are mostly empty of energy. This simulated signal is nice for illustrating how spectral decompositions works but is not representative of realistic scenarios. Most real signals involve high-frequency components riding on top of low frequency signals, so let’s explore another simulated signal with a temporally limited burst of high-frequency “noise” added to a low-frequency signal. The code for simulating the signal and plotting the spectrogram and scalogram is as follows:   proc iml; /*We start by simulating a digital signal with a short burst of high-frequency noise*/ f_s = 500; /*sampling frequency*/ f_s2 = f_s/2; /*Nyquist frequency*/ t_s = 1/f_s; /*sampling period (time between samples)*/ N_samp = 1000; /*number of samples*/ t1 = do(0, (N_samp/2-1)*t_s, t_s); t2 = do((N_samp/2)*t_s, (N_samp/2+5)*t_s, t_s); t3 = do((N_samp/2+6)*t_s, (N_samp-1)*t_s,t_s); f1 = 10; f2 = 200; pi = constant("pi"); x1 = sin(2*pi*f1*t1); /*signal for the first half*/ x2 = 0.1*sin(2*pi*f2*t2) + sin(2*pi*f1*t2); /*short piece of HF noise*/ x3 = sin(2*pi*f1*t3); /*signal for second half*/ t = t1||t2||t3; /*combined time from N=0 to N=200*/ x_in = x1||x2||x3; /*combined signal with short HF noise*/ create work.input_signal var {"t","x_in"}; append; close work.input_signal; /*plot our signal in the time domain*/ submit; title "Input Digital Signal"; proc sgplot data=work.input_signal; series x=t y=x_in; xaxis grid label="Time (s)"; yaxis grid label="Amplitude"; run; endsubmit; /*Now we obtain our signal's frequency domain representation using STFT*/ winLen = 50; /*this will have a big impact on the time/frequency resolution*/ window = tfwindow(winLen, "Gaussian"); noverlap = 0; nfft = 64; /*increasing this increases the number of frequency bins*/ ramp = palette("Spectral",7)[7:1]; /* reverse "spectral" color ramp */ title "STFT Spectrogram"; call spectrogram(x_in,window,noverlap,nfft) scale="log" colorramp=ramp sampRate=f_s panel=1; x_in = x_in`; /*initialize the wavelet graphics module*/ %wavinit; /*set wavelet options using macro variables from %wavinit macro*/ /*in this case let's use a predefined Daubechies 3-tap wavelet*/ optn = &daubechies3; /*now perform the discrete wavelet transform to find the decomposition*/ call wavft(decomp, x_in, optn); /*plot a scalogram of the wavelet decomposition*/ title "Wavelet Transform Scalogram"; call scalogram(decomp); quit;   This yields the following spectrogram and scalogram:     We can see in the spectrogram that there is some amount of high-frequency noise added to the signal at about 1 second. The noise is low power, so it’s exact frequency is ambiguous in the spectrogram. The width of the time window also makes it hard to pinpoint when this noise was added to the signal.     In the wavelet scalogram we still don’t have good frequency resolution at the higher frequencies, but we can see that there was some high-frequency noise added about halfway through the signal. At the 9 th level we can see a very narrow yellow line, corresponding to the exact portion of the signal where the high-frequency noise was added. The wavelet scalogram does a much better job of localizing the high-frequency noise in time than the STFT spectrogram.   Wavelets are a powerful tool for analyzing signal and sensor data. We have seen how we can use wavelets to decompose periodic signals into time-frequency representations that often have more useful resolution properties than the STFT spectrogram, especially for real-world signals. The examples in this blog just scratch the surface of what we can do with wavelets, and we will explore applications of wavelets such as denoising and data compression (for one-dimensional signals and two-dimensional images) in future blog posts.   References: Previous article on time-frequency analysis using Short-Time Fourier Transform SAS Documentation describing wavelet analysis SAS Documentation page listing the wavelet options and symbolic names   ... View more

The purpose of this post is to learn how to detect anomalous frequencies in streaming signal data using SAS Event Stream Processing. We learn how to implement the streaming Short-Time Fourier Transform algorithm using the Calculate Window, and we also learn how to explore and interpret the results of the streaming algorithm. Illustrative screenshots of the Event Stream Processing Studio graphical interface are included throughout the post, and a full XML Project containing all the work described in this post is provided at the end as a reference.   Deploying Models in SAS Event Stream Processing   SAS Event Stream Processing allows us to analyze real-time data streams and immediately react to patterns or information in the data stream. This is useful when working on anomaly detection problems because we usually want to take some action after finding anomalies, and our time to take this action can be limited. If we detect an anomaly in the manufacturing line, we may want to shut down the line before damage to the products starts to accumulate. If we react too late to the anomaly we may lose a lot of our production to damage, so the faster we can identify the anomaly and shut down the manufacturing line, the more profitable our business will be. SAS Event Stream Processing provides the ability to deploy an anomaly detection model in real-time, so we can go immediately from streaming in data to detecting anomalies and reacting to them. We can even use SAS Event Stream Processing to send a signal to the manufacturing line to shut down immediately upon detection of an anomaly (assuming we can interact with the manufacturing line via an API). Previous posts on the topic of anomaly detection provide some background for the concepts and algorithms discussed in this post. In particular, the posts linked below introduce some basic approaches to anomaly detection and the Short-Time Fourier Transform algorithm:   Introduction to Supervised and Unsupervised Anomaly Detection Identifying Anomalous Frequencies in Signal Data with Short-Time Fourier Transform using SAS/IML   Here we will be focused on using SAS Event Stream Processing software to implement the Short-Time Fourier Transform on streaming signal data.   Real-Time Detection of Anomalies using Short-Time Fourier Transforms   The short-time Fourier transform allows us to decompose signals into frequency components over time. This decomposition can reveal when unwanted frequency components appear in the data. A common example of this is the unwanted appearance of high-frequency noise in the data. This could be the result of unwanted vibrations in mechanical devices, or noise coupling in circuits. Regardless of how the noise enters the system, we want to identify when the noise appears and flag the times where the signal contains the unwanted noise. We may also want to use digital filters to remove the unwanted noise, but for now we will focus on detecting the anomalies.   In a previous post (Identifying Anomalous Frequencies in Signal Data with Short-Time Fourier Transform using SAS/IML) we learned about using the STFT in SAS/IML using the SPECTROGRAM CALL method. Here we will follow up on that problem and deploy our anomaly detection method in real-time using SAS Event Stream Processing. As a reminder we had a 30 Hz signal with high-frequency noise added for a very short time, and we want to detect when this high frequency noise appears. The spectrogram we plotted for this signal revealed the high-frequency noise was localized around 200 Hz right after the 1 second mark in the signal:   Select any image to see a larger version. Mobile users: If you do not see this image, scroll to the bottom of the page and select the "Full" version of this post.   The patch of orange points around 200 Hz after 1.0 seconds reveals the anomalous frequency added to the data. We can characterize the anomaly as any frequency component in the data above about 100 Hz (accounting for spectral leakage from the normal 30 Hz signal) with non-negligible power. Looking at the power spectrum it looks like above 100 Hz the log(Power) values are all below -4 except for the anomalous frequency components. We can choose a threshold value for our anomaly detection model, identifying any frequency component above 100 Hz with log(Power) greater than -4 as an anomaly. The basic idea is that the 30 Hz signal is the normal data, and anything that deviates significantly from the normal data (like the burst of 200 Hz noise we see in the spectrogram) is an anomaly. Our goal is to output real-time information about when the anomaly occurs in the data and at which frequencies it occurs (knowing that in this example frequencies below 100 Hz are non-anomalous).   We start by creating a project in SAS Event Stream Processing Studio and adding some relevant windows to the project. As we walk through the project, we will see screenshots of portions of the project in SAS Event Stream Processing Studio. This is included as a reference to help readers learn to use SAS Event Stream Processing Studio, but it is not provided as a method to reproduce the steps in this post. SAS Event Stream Processing Studio projects are based on XML, so behind the scenes in the graphical interface, XML is created to represent the project. We have provided the full XML for each of the examples discussed in the post at the end in the Reference Materials section. To explore these projects on your own, you can import the XML into SAS Event Stream Processing Studio and navigate the project using the graphical interface.     The first window in the project (readDigitalSignal) is Source Window and is used to read the input data. In real project this would read data from a sensor device or a website that is designed to output streaming data, but in this toy example we read the signal data from a CSV file. We use the File/Socket Connector to read the CSV file into the source window, and the source window will output the information read from the CSV file, which in this case is just the time variable t, and the input signal, x_in. The output schema for the source window shows that the time variable t is used as the Key Field for the streaming project.     The second window in the project (shortTimeFourierTransform) is a Calculate Window and is used to perform the STFT on the streaming data. The calculate window will use the input streaming data (in this case the time variable t and the input signal, x_in) and the STFT algorithm to calculate the power across frequency bins in each time window. The STFT settings like the window length/type and the overlap between windows can be set using the SAS Event Stream Processing Studio graphical interface (hardcoded), or they can be set by referring to an input configuration file living on the operating system. In this example we hardcode the STFT parameter settings:     The advantage of using a configuration file is the ability to modify settings without altering the project, allowing the same anomaly detection project to be reused in different settings and environments. The calculate window will output the power and phase in each frequency bin for each time window. Notice in the output schema that the combination of the time and bin fields is used as the Key Field:     A key requirement for using the STFT algorithm in SAS Event Stream Processing is to specify a map between the STFT algorithm and the input/output fields in the schema. The input map is simple, we just specify that the field t is used as the time ID variable and the field x_in is used as the input signal:     For the output map we must specify where we store the STFT results, but this time it can be a bit confusing since there are two options for how we store the outputs, a list format output, and a non-list format output. We use the non-list format output in this example:     In the output map above, we leave the keyOut Role blank because we use a combination of the timeIdOut and the binOut variables as the key variable (as we saw in the output schema earlier). We also leave the powerListOut and phaseListOut Roles blank because we are not using the list format. This is important since we will get a confusing error if we try to run the pipeline with all the output roles specified. Using the Output Map above we generate an output event for each time/bin combination, meaning with 32 bins we will get 32 output events for each time point. The list format instead outputs a single event for each time point, with a list of power and phase values for each frequency bin. This is a minor difference in format, but when configuring the SAS Event Stream Processing project, it is important to choose a single output format and stick with it to avoid errors and confusion.   The third window in the project (binToFreq) is a Compute Window used to convert the numeric frequency bins (just an integer number counting the bins from 1 to 32 in this example) into actual frequencies in Hertz. We do this by multiplying the bin number by the sampling frequency and then dividing by the total number of bins in the schema. This creates a new field in the data containing the frequencies in Hertz.   The fourth and final window in the project (noiseDetection) is just a Filter Window that filters the incoming stream to only output events when anomalous signals are detected, in this case it detects the presence of high frequency noise (signal components above 100 Hz with power greater than 0.01). The output schema for the filter window is the same as the input schema, we just filter out events based on the expression specified below in the filter window (notice how we use the field freq created by the compute window as a filter variable):     In this simple example we just filter the noise out from the signal using rules we learned from examining the spectrogram, but we could also use more sophisticated tools for finding peaks in the spectrogram to identify anomalous frequencies. SAS Event Stream Processing has a peak finding algorithm in the Calculate Window that can be used to look for peaks in plots like the spectrogram and then output the time at which these peaks occurred and the signal amplitude at the peak. In a real deployment we would want to output this detected anomaly to another downstream system, either by sending an API call or message to the downstream system, or by writing output data (like a CSV file) that is consumed by a downstream system.   We can test the SAS Event Stream Processing project by saving it and then selecting Enter Test Mode in SAS Event Stream Processing Studio. This allows us to see how the anomaly detection pipeline will work using the simulated signal data, and it will allow us to test our project to make sure all the settings work. Once in test mode we can run a test deployment to see how the CSV data looks when it is streamed through the pipeline:     We can see that the source window reads in the data from the CSV file and assigns an Opcode to each observation in the data, creating events in SAS Event Stream Processing.     The Calculate window outputs the STFT results, which in this case is an event for each combination of time and frequency bin (if we were using the list format for the STFT algorithm we would get different results). Notice that the power in each frequency bin is output in a linear scale, normally when we plot spectrogram results, we use a log or decibel scale for spectral power. This is something to keep in mind when comparing the output of the STFT algorithm in SAS Event Stream Processing to the output we would see when using the call spectrogram routine in SAS/IML.     The binToFreq Compute Window just converts the bin number to a frequency in Hz, and the noiseDetection Filter Window is used to select only those events that contain unwanted high-frequency noise. We can see that while the previous windows had 1000 rows each, the noiseDetection window only has 22 rows, corresponding to the 22 output events where high frequency noise was detected. This doesn’t mean 22 time points contained high frequency noise, since each time point can have multiple events corresponding to the different frequency bins with unwanted noise. At this point in the project, we could use a subscriber connector to write out the results or to connect to another downstream system that uses the time information in the output to perform some kind of business task. A simple example would be a system where the input signal is voltage measured on a circuit in a wind turbine, and when high frequency voltage is measured in the circuit, we send a signal to the circuit to connect to a ground and discharge the high-frequency voltage. In this example we are trying to protect the turbine generator from the negative effects of unwanted voltage components while still maximizing the uptime of the turbine power generation.   This is just one way we can apply the STFT algorithm to input signal data (we can also use SAS/IML and SAS Visual Forecasting procedures), but it is the best way to perform a STFT on streaming data, especially if we need to react to the anomalies in a time-sensitive way. This demonstration does not include the final step of reacting to the anomalies in real-time, but in a situation where we are detecting anomalies in a manufacturing setting, it is usually as simple as sending a shutdown signal to a programmable logic controller (PLC) device on the manufacturing line. This is also just one algorithm we can use to detect anomalies in real-time, a subsequent post will discuss deploying support vector data description models in SAS Event Stream Processing. This differs from the current approach because we will need to deploy a model using a scoring artifact (an ASTORE file) instead of just using an online algorithm in the Calculate Window.   Reference Materials - STFT_Anomaly_Detection SAS Event Stream Processing Project XML:   <project name="STFT_Anomaly_Detection" threads="1" pubsub="auto" heartbeat-interval="1"> <metadata> <meta id="studioUploadedBy">student</meta> <meta id="studioUploaded">1685034102948</meta> <meta id="studioModifiedBy">student</meta> <meta id="studioModified">1686236532700</meta> <meta id="layout">{"cq1":{"binToFreq":{"x":-80,"y":20},"noiseDetection":{"x":-80,"y":155},"readDigitalSignal":{"x":-80,"y":-255},"shortTimeFourierTransform":{"x":-80,"y":-120}}}</meta> </metadata> <contqueries> <contquery name="cq1"> <windows> <window-source index="pi_EMPTY" insert-only="true" name="readDigitalSignal"> <description><![CDATA[This window reads in the digital signal for analysis in SAS Event Stream Processing. This toy example uses a CSV file, but we will illustrate reading data from a URL as well.]]></description> <schema> <fields> <field name="t" type="int64" key="true"/> <field name="x_in" type="double"/> </fields> </schema> <connectors> <connector class="fs" name="Signal_In_CSV"> <properties> <property name="type"><![CDATA[pub]]></property> <property name="header"><![CDATA[1]]></property> <property name="addcsvopcode"><![CDATA[true]]></property> <property name="addcsvflags"><![CDATA[normal]]></property> <property name="fsname"><![CDATA[/mnt/data/anomaly/simulated_signal_data_w_noise.csv]]></property> <property name="fstype"><![CDATA[csv]]></property> </properties> </connector> </connectors> </window-source> <window-calculate algorithm="STFT" index="pi_EMPTY" name="shortTimeFourierTransform"> <description><![CDATA[This window takes the input digital signal and performs a short-time Fourier Transform to identify the frequencies in the signal as a function of time.]]></description> <schema> <fields> <field name="time" type="int64" key="true"/> <field name="bin" type="int64" key="true"/> <field name="power" type="double"/> <field name="phase" type="double"/> </fields> </schema> <parameters> <properties> <property name="windowLength"><![CDATA[25]]></property> <property name="windowType"><![CDATA[15]]></property> <property name="windowParam"><![CDATA[-1.0]]></property> <property name="fftLength"><![CDATA[64]]></property> <property name="overlap"><![CDATA[10]]></property> <property name="binsInSchema"><![CDATA[32]]></property> </properties> </parameters> <input-map> <properties> <property name="input"><![CDATA[x_in]]></property> <property name="timeId"><![CDATA[t]]></property> </properties> </input-map> <output-map> <properties> <property name="timeIdOut"><![CDATA[time]]></property> <property name="binOut"><![CDATA[bin]]></property> <property name="powerOut"><![CDATA[power]]></property> <property name="phaseOut"><![CDATA[phase]]></property> </properties> </output-map> <connectors> <connector class="fs" name="STFT_Write_CSV"> <properties> <property name="type"><![CDATA[sub]]></property> <property name="header"><![CDATA[full]]></property> <property name="snapshot"><![CDATA[false]]></property> <property name="fsname"><![CDATA[/mnt/data/anomaly/simulated_signal_STFT.csv]]></property> <property name="fstype"><![CDATA[csv]]></property> </properties> </connector> </connectors> </window-calculate> <window-filter index="pi_EMPTY" name="noiseDetection"> <description><![CDATA[This filter window looks for non-negligible power in the high frequency bins. Spectral leakage means power is spread out across many bins so we look for non-negligible power (above 0.01 in the linear scale) in frequency bins above 100 Hz.]]></description> <expression><![CDATA[(freq > 100) and (power > 0.01)]]></expression> </window-filter> <window-compute index="pi_EMPTY" name="binToFreq"> <description><![CDATA[The STFT outputs information about the power in different frequency bins. The actual frequency corresponding to each of these bins is based on the number of bins and the sampling frequency. frequency = (bin / binsInSchema) * nyquist_frequency]]></description> <schema> <fields> <field name="time" type="int64" key="true"/> <field name="bin" type="int64" key="true"/> <field name="power" type="double"/> <field name="phase" type="double"/> <field name="freq" type="double"/> </fields> </schema> <output> <field-expr><![CDATA[power]]></field-expr> <field-expr><![CDATA[phase]]></field-expr> <field-expr><![CDATA[(bin * 250) / 32]]></field-expr> </output> </window-compute> </windows> <edges> <edge source="readDigitalSignal" target="shortTimeFourierTransform" role="data"/> <edge source="shortTimeFourierTransform" target="binToFreq"/> <edge source="binToFreq" target="noiseDetection"/> </edges> </contquery> </contqueries> </project>     Find more articles from SAS Global Enablement and Learning here. ... View more

The purpose of this post is to learn how digital signal processing tools can be used to identify anomalous frequencies in sensor data. We start by learning about digital signals and the Fourier Transform before introducing the short-time Fourier Transform (STFT). We then use the SAS Interactive Matrix Language (SAS/IML) to simulate a digital signal with unwanted high-frequency noise and use the STFT to identify the presence of this unwanted noise.   Introduction to Analog and Digital Signals   Sensors can provide useful real-time information about the condition or quality of devices or components. These sensors could measure simple scalar properties like voltage in a circuit or temperature in a kiln, but they could also measure more complex vector quantities like the wind speed and direction at a manufacturing facility. In both cases the sensor will output data that effectively looks and acts like time series data (each observation output by the sensor will have a value for when it is output and a scalar or vector value for the measure of interest), but it is more useful to think of this sensor data as a digital signal instead of treating it as a time series.   Digital signals are sampled versions of continuous analog signals. Analog signals are continuous functions of time and cannot be represented in a computer. A very simple example of an analog signal would be the voltage output from an ideal alternating current (AC) outlet as measured by an analog oscilloscope:   Select any image to see a larger version. Mobile users: If you do not see this image, scroll to the bottom of the page and select the "Full" version of this post.   Some important properties of periodic signals like this are the amplitude of the signal, the frequency of the signal, and the phase of the signal. The amplitude defines the height of the peaks in the signal, while the frequency is the inverse of the time between the peaks (the time between the peaks is called the period of the signal). Of course, most real signals are not just simple sine waves like the ideal AC outlet and often contain multiple components with different frequencies and amplitudes.   Although a sensor device may measure the voltage in a circuit every 0.5 seconds and output a series of numbers corresponding to the voltage over time, the actual voltage over time is a continuous function of time instead of a collection of discrete samples. The process of sampling the circuit every 0.5 seconds turns the analog signal into a digital signal; if this process is done carefully the digital signal will represent the phenomenon of interest in the analog signal. The key idea is that the sampling frequency (i.e., the frequency with which we take digital samples from the underlying analog signal) sets an upper limit on the frequencies in the original analog signal that can be identified in the sampled digital signal. The Nyquist frequency is equal to half the sampling frequency and is the largest frequency that can be represented correctly by the digital signal. This means when sampling analog signals, we must choose the sampling frequency to be at least twice as large as the largest frequency in the signal, else we will lose important information moving from the analog signal to the digital signal.     In picture above the orange dots are the samples in the digital signal while the blue line shows the original analog signal. A good sampling frequency like the one above successfully captures a representation of the analog signal with the samples.     Now in the picture above we try to use a lower sampling frequency with the same signal. This time our equally spaced sampling points always correspond to zeros in the analog signal. The digital signal sampled here will just be a collection of samples all with an amplitude of zero, even though the analog signal was not just constantly zero. A bad sampling frequency yields unclear results and fails to capture the information present in the original analog signal.   Introduction to Fourier Transforms   Since most real-world signals contain multiple components with different frequencies and amplitudes, just looking at the numeric values of the digital signal won’t provide much insight or understanding about the device that generated the signal. It is often more useful to break down the signal into its different frequency components and then analyze those components to either gain an understanding about the device, or to create inputs for predictive models. This conversion from the time domain to the frequency domain can be accomplish using the Fourier Transform, a mathematical approach to turn a signal as a function of time into a signal as a function of frequency. When working with digital signals we use a Discrete Fourier Transform (DFT) since we are operating on discrete samples, and it turns out there is an incredibly efficient algorithm for calculating DFTs called the Fast Fourier Transform (FFT). The FFT is one of the most important algorithms of the past century and enabled a lot of the functionality we are used to with modern digital computers and communication devices. The fundamental idea behind the Fourier transform is best illustrated in pictures:     This is a plot of digital signal in the time domain: the x-axis is time in seconds, and the y-axis is amplitude. This signal only contains a single frequency component, and it looks like the time between peaks is about 0.2 seconds, so the frequency would be (0.2) -1 = 5 Hz. We can see this by looking at a plot of the Fourier transform of this signal (calculated using the FFT algorithm):       This is a plot of the digital signal in the frequency domain: the x-axis is frequency in Hertz (inverse seconds) and the y-axis is amplitude. Although the peak is not a perfect vertical line (spectral leakage and limited frequency resolution leads to non-ideal results, but we won’t get too detailed in this blog post), we can see that there is a peak at 5 Hz, matching what we expected from looking at the signal in the time domain. This simple example illustrates the basics of the output we get from Fourier transforms, but we don’t gain much from looking at this simple signal in the frequency domain. The story is changed when we have multiple different frequency components in the signal.     In this example our signal has two different frequency components, with short peaks mixed in with taller peaks in the time domain. Knowing our signal has two frequency components we can probably estimate the different frequencies by looking at the time between the peaks, but it is much easier to interpret the signal contents when we take a Fourier transform and move to the frequency domain:     Here we can clearly see that there is a small peak at 5 Hz, and a larger peak at 10 Hz. Real-world signals will contain many different frequency components and so the Fourier transform of these signals will have many peaks in the frequency domain, some of which may even overlap.   Introduction to Time-Frequency Analysis and the Short-Time Fourier Transform   The Fourier transform makes it easy to analyze the components of signals with constant frequency, but if the frequency changes during the signal, the Fourier transform will not tell us when this change occurs (or even that it occurs at all), instead it will just output all the frequencies in the signal at any time in the signal duration. Basically, the Fourier transform is time-independent; if we want to understand the time-dependence of frequencies in the signal we must use a different approach. This is best illustrated by looking at the Fourier transforms for two very different signals with the same basic frequency components:     This first signal starts with a 5 Hz sinusoid and after 1 second the frequency increases to 15 Hz. When we take the Fourier transform of the signal using the FFT algorithm, we see two peaks in the frequency domain, one for the 5 Hz portion of the signal and one for the 15 Hz portion of the signal. Although this tells that the signal has two different frequency components, it gives us no information about the timing, so we are not able to see that the signal only has one frequency at a time, it just changes after 1 second.     This second signal is just a linear combination of a 5 Hz sinusoid and a 15 Hz sinusoid. There is no time dependence in the signal and the signal always has these two frequency components. The Fourier transform basically shows the same things as the previous example, a peak at 5 Hz and a peak at 15 Hz. In these examples we have simple signals, so we can look at the time domain to distinguish between the signal with the time-dependent frequency and the signal with the time-independent frequencies, but the frequency domains for these two signals look the same. When working with more complex real-world signals we won’t find it easy to interpret frequency changes in the time domain, and using the Fourier transform to move to the frequency domain won’t help in identifying timing information.   The solution to our problem is to perform a sequence of windowed Fourier transforms, with the windows each covering a small portion of time, so that we can plot a sequence of Fourier transforms over time to uncover how the frequencies change. This approach is called a Short-Time Fourier transform (STFT) and the result is a collection of individual Fourier transforms for different time windows of the original signal.     Instead of performing a Fourier transform on the entire signal, we break it into two different windows, performing a Fourier transform on the first window from 0-1 second and then a performing a second Fourier transform on the window from 1-2 seconds. Each of these transforms will yield a frequency spectrum, which can then be plotted in sequence. This approach reveals that the signal only contains a 5 Hz component from 0-1 second and only contains a 15 Hz component from 1-2 seconds.   More realistic STFT calculations involve many more windows and are often plotted using a two-dimensional heatmap, but the fundamental idea is the same as in this toy example. We can access FFT and STFT functionality in a couple of different ways using SAS software. The SAS Interactive Matrix Language (SAS/IML) contains a variety of easy-to-use digital signal processing tools, including a routine to calculate FFTs and STFTs. The TSMODEL procedure in SAS Visual Forecasting (SAS VF) also includes a Time-Frequency Analysis (TFA) package with objects for performing FFTs and STFTs. Finally, SAS Event Stream Processing (SAS/ESP) includes the STFT as an online algorithm that can be applied to windowed streaming data. Let’s explore an example of plotting the results of a STFT using SAS/IML. In this context SAS/IML is most useful for exploring digital signals and prototyping workflows, SAS VF is most useful for analyzing traditional time series using digital signal processing tools, and SAS/ESP is most useful for deploying digital signal processing workflows.   Spectrogram STFT Example using SAS/IML   To illustrate the use of STFT for anomaly detection we simulate a simple sinusoidal signal with a short period of injected low-amplitude, high-frequency noise in the middle of the signal duration. To be specific, we simulate a 30 Hz signal for 1000 time points (with a sampling frequency of 500 Hz this corresponds to a signal with a 2 second duration), and we inject 50 samples of 200 Hz noise for a duration of 50 samples (0.1 seconds). This noise is injected with an amplitude of 5% the original signal amplitude, so it doesn’t make a noticeable change to the time-domain plot of the signal, but its presence becomes noticeable in the time-frequency domain. This example is generic by intention to illustrate the core concepts of STFT, but we could imagine this signal being a voltage measurement from a circuit that has temporarily been coupled to a source of high-frequency noise, or a sound/vibration measurement from a turbine device with a damaged blade that scrapes against a surface briefly during the operation. Either way, our example is a situation where we ideally expect a consistent 30 Hz signal to be output from our device, but some anomalous event has injected unwanted high-frequency noise in our signal, and we want to identify when this happens so we can potentially remove it.   We start in SAS/IML by defining the sampling frequency and establishing an array of time values for the simulated signal. The detailed code below to create the data is included as a reference for interested readers (and to make it easier to reproduce the demonstration). If you are not as comfortable with SAS/IML syntax, the main focus of this blog is on using the Spectrogram subroutine to visualize the frequency spectrum. The code to plot the spectrogram is the last block of code after the plot of the input signal, and although it is still part of SAS/IML, the syntax is fairly simple.   proc iml; f_s = 500; /*sampling frequency*/ f_s2 = f_s/2; /*Nyquist frequency*/ t_s = 1/f_s; /*sampling period (time between samples)*/ N_samp = 1000; /*number of samples*/ /*array of times (first half)*/ t1 = do(0, (N_samp/2-1)*t_s, t_s); /*array of times (section with high-frequency noise)*/ t2 = do((N_samp/2)*t_s, (N_samp/2+50)*t_s, t_s); /*array of times (second half)*/ t3 = do((N_samp/2+51)*t_s, (N_samp-1)*t_s, t_s);   The sampling frequency is 500 Hz, so the largest frequency we can identify in the digital signal will be 250 Hz. The arrays t1 and t3 will be used to simulate the normal signal, while the 50-sample array t2 will be used to simulate the signal with the high-frequency noise injected.     f1 = 30; f2 = 200; pi = constant("pi"); /*signal for the first half*/ x1 = sin(2*pi*f1*t1); /*signal corrupted by a burst of high-frequency noise*/ x2 = sin(2*pi*f1*t2) + 0.05*sin(2*pi*f2*t2); /*signal for the second half*/ x3 = sin(2*pi*f1*t3); /*combined time from N=0 to N=2000*/ t = t1||t2||t3; t = t*1000; /*convert time to milliseconds*/ /*combined signal with a burst of high-frequency noise*/ x_in = x1||x2||x3; create work.input_signal var {"t","x_in"}; append; close work.input_signal;   The signals x1 and x3 are the normal signal with frequency 30 Hz, while the signal x2 is the short duration with the 5% amplitude high-frequency noise injected into the signal. We can see that x2 is just a linear combination of the original normal signal with the low-amplitude noise. After defining the 3 different time arrays and 3 different signal arrays, we concatenate them in SAS/IML to create the overall time array t and the overall signal array x_in. We write these arrays to a SAS dataset called work.input_signal so we can plot them using SGPLOT.   /*plot our signal in the time domain*/ submit; title "Input Digital Signal"; proc sgplot data=work.input_signal; series x=t y=x_in; xaxis grid label="Time (ms)"; yaxis grid label="Amplitude"; run; endsubmit;     Plotting the signal as a function of time we can see there looks like there is some distortion around 1000 ms in the signal, but it’s unclear exactly what is happening in the signal, and it doesn’t look like there is a 200 Hz component added into the signal. Visually, the frequency of the signal in the above plot looks mostly constant; of course, visually inspecting the plot of a signal is not a good way to determine the frequency components in the signal. Let’s use the STFT to plot a spectrogram.   /*The window length will have a big impact on the time/frequency resolution*/ winLen = 25; window = tfwindow(winLen, "Hanning"); noverlap = 10; /*increasing nfft zero-pads the windows and increases the number of frequency bins*/ nfft = 64; /* reverse "spectral" color ramp */ ramp = palette("Spectral",7)[7:1]; title "STFT"; call spectrogram(x_in,window,noverlap,nfft) scale="log" colorramp=ramp sampRate=f_s panel=1; quit;   Before we look at the spectrogram output lets briefly discuss some of the important STFT options used in this example:     Window Length – the window length (winLen = 25) determines how many time samples are included in each time window of the STFT. Longer window lengths mean more samples per window and thus higher frequency resolution. Shorter window lengths mean more windows overall given a fixed signal duration and thus higher time resolution. There is always a fundamental trade-off in resolution between time and frequency when performing the STFT. Window Type – the window function choice affects spectral leakage across the frequency bins. In this example we use a Hann (often called Hanning) window function for each segment. The basic idea with windowed transforms is to eliminate negative effects from using a rectangular window (which is used by default if another window function is not chosen). The details of how and why this works are outside the scope of this blog, but spectral leakage and window functions are an important topic in digital signal processing. FFT Length – The number of FFT samples (nfft = 64) is the number of points used by the FFT algorithm to calculate the frequency spectrum for each time window in the STFT. This number must be greater than or equal to the window length; if nfft = winLen then no zero-padding will occur and the FFT will be performed on the windowed signal. If nfft > winLen then the windowed signal will be padded with zeros until it has the length specified by nfft. These added zeros help by ensuring we can use efficient FFT algorithm (by ensuring that the number of samples is a power of 2), but they also help by increasing the number of frequency bins in the FFT output. This doesn’t fundamentally increase the frequency resolution of the STFT since the output will be interpolated across the frequency bins, but it does improve the plotting of the results. Window Overlap – the overall between windows (noverlap = 10) determines how many samples overlap between subsequent windows in the STFT analysis. Choosing no overlap would mean that information near the boundaries of each window is lost, while choosing a very large overlap would mean a lot of redundant information is included in each window, increasing processing time unnecessarily. In general, a good practice is to choose the overlap to be around 50% of the window size. In this example we use an overlap of 10 with a window size of 25, so our overlap is 40% of the window size. We use the call spectrogram() routine in SAS/IML to calculate the STFT and plot it in a two-dimensional heatmap, all in one function call. We provide the input signal x_in, the STFT settings described above, the sampling frequency for the underlying signal so we can output frequency in Hertz, and a scale for the power in each frequency bin. We use the log scale (the decibel scale is also useful here) since the power values outside of the 30 Hz signal are hard to distinguish in the linear scale. The option “panel=1” is chosen so that we output a panel of results, including the time domain representation of the signal, the time-frequency heatmap, and an averaged view of the frequency spectrum across all time.     Looking at the output of the spectrogram, we can see the 30 Hz component contains most the power in the heatmap. On a linear scale this would dominate all the other frequency components, making it hard to see lower power high-frequency components in the plot. On the log scale we can see there is non-negligible power around the 150-250 Hz range between Time = 1.0 seconds and Time = 1.1 seconds. The orange smear around 200 Hz corresponds to the injection of the high-frequency noise into the signal, and the spectrogram provides us detailed information about the timing of when the noise appeared in the signal, along with the frequency of the injected noise. Fundamental limitations associated with the STFT still apply, so we don’t perfectly localize the noise in time and frequency. It looks like the high-frequency noise is spread out between 175 and 225 Hz (due to spectral leakage), and the pixelation of the spectrogram clearly illustrates the limited resolution in time and frequency. We could improve the frequency resolution by using longer windows (and thus fewer windows) in the STFT, but this would degrade the time resolution, making it more difficult to identify when the noise was injected. Here is an example of changing some settings as follows:   winLen -> 250 noverlap -> 100 nfft -> 256     Here we can see that the frequency bands are much narrower, although we still have some spread around the true frequencies of 30 Hz and 200 Hz. The tradeoff is in the time resolution, it looks like the high-frequency noise was added around 0.75 seconds and didn’t end until about 1.25 seconds. This loss in time resolution makes sense since we used longer windows, so the high-frequency noise appears (with at least a little bit of power) in all the time windows between 0.75 and 1.25 seconds.   Anomaly Detection Using STFT   We can use the STFT to perform anomaly detection by collecting normal signals (signals that we are confident do not contain anomalies) and looking at their frequency content. We plot the spectrogram of the non-anomalous data and identify major frequency components and when they occur in the signal. This characterizes the normal data, which we will use as a baseline for judging whether signals contain anomalous frequencies. We could then perform a real-time STFT on all new signals and compare the frequencies we see in the spectrogram with those from the normal FFT. Any extra frequency components we identify in the new signals would correspond to anomalies in the data (assuming we did an adequate job of characterizing the normal data signals). This is a kind of unsupervised anomaly detection where we don’t have to collect data representing the anomalies, we only collect data to represent the normal operating conditions of the system that generates the signal.   Continuing with the example of a 30 Hz signal with 200 Hz noise injected, our normal signal would be just the 30 Hz signal with no 200 Hz noise.     The spectrogram of the normal operating conditions reveals that most of the power is concentrated around 30 Hz. To be more specific, we can identify that the normal operating conditions have no significant frequency components above 100 Hz (the 30 Hz signal leaks spectral information around the true frequency), so anytime we see significant power above 100 Hz we should identify it as an anomaly. Of course, we must be specific about what significant power means, so we are looking for signals with the log of the power greater than about -4 in frequency ranges above 100 Hz. Since we want the log of power to be greater than -4, on the linear scale we want power > 0.018; this can be used to threshold anomalies when we go to evaluate new signals for anomalies using the spectrogram.   This approach for detecting anomalies can be done entirely using SAS/IML using the digital signal processing tools, but we could also use the SAS Visual Forecasting tools in the Time Frequency Analysis Package (TFA) in PROC TSMODEL to perform the STFT to detect anomalies. SAS/IML and the TFA package are great for offline analysis of digital signals, but to perform real-time scoring of anomalies we would need to use the Calculate Window with the STFT algorithm SAS Event Stream Processing. A subsequent blog post will discuss deploying anomaly detection models in SAS/ESP.   References/Useful Links:   Time-Frequency Analysis using SAS® Getting Started with the SAS/IML® Language CALL SPECTROGRAM SAS/IML SAS Help Center: Digital Signal Processing The Fourier Transform Scientific American Article by Ronald N. Bracewell This is a great introductory article talking about Fourier analysis and the FFT algorithm     Find more articles from SAS Global Enablement and Learning here. ... View more

The purpose of this blog is to explain the basic concepts of anomaly detection and to compare supervised and unsupervised anomaly detection methods. This blog will use a supervised gradient boosting model (PROC GRADBOOST) and an unsupervised support vector data description model (PROC SVDD) to detect anomalies in a sample dataset.   Basic Introduction to Anomaly Detection   Anomalies are rare events that deviate significantly from most of the data. These anomalies could be extreme events that are generated in the same ways as the normal data (extreme weather events are caused by the same underlying process as any other weather event), but they could also be events that are generated in a completely different way from the normal data (fraudulent bank transactions are not performed for the same reasons as normal bank transactions). Detecting these anomalies in data can be useful for a variety of reasons, and the value in detecting them will be different depending on the type of data: Anomalies in data can indicate that our model for the underlying data generation process is inadequate or incorrect. Example: Rogue waves in the ocean are unusually high waves and are extremely unlikely when representing the ocean surface as Gaussian field. This is because the Gaussian Sea model for ocean sea surface is based on linear waves and rogue waves are a result of nonlinear effects in the Navier-Stokes equations. Anomalies in data can represent unwanted or malicious activity that must be addressed in a specific way. Example: Banks lose money on fraudulent transactions, so identifying fraud in advance and preventing these transactions is a profitable activity for banks. Of course, fraud is profitable for the fraudsters as well so they will continually change their behavior, meaning the anomalies in the transaction data will change over time as the fraudster adapt to the bank’s fraud prevention. Anomalies in machine sensor data can provide a leading indicator for component failure. Example: A rotating blade used to cut plastic on a manufacturing line outputs torque and temperature information each second. At some point during the operation the torque and temperature measurements start to increase beyond the normal operating range. This anomaly is a leading indicator that the blade is scraping against a surface in the manufacturing line, leading to larger torque to get past the surface and increased temperature due to friction. Turning off the manufacturing line early and replacing the blade before it breaks and ruins the production can save time and money. Anomalies can represent outliers that should not be included in data used for predictive models or statistical analysis. Example: A grocery store offers a rewards card account and collects data about the purchase history for each customer account and plans to use this data for targeted marketing. Some of the accounts have extremely large monthly spending totals, indicating that these customers are big spenders at the store (and that marketing to them might be effective). A more detailed look at the data reveals that these large monthly spending totals are associated with accounts used by cashiers when customer forget to bring their rewards card or cannot remember their account number. These “big spenders” would not be useful to include in an analysis of customer spending habits since they don’t represent a customer. One key takeaway from these different examples is that the goal in anomaly detection isn’t always to just remove the anomalies from the data, sometimes we are more interested in the anomalies than in the normal data. This contrasts with the traditional statistical treatment of outliers, which is usually to exclude them from the data analysis. Sometimes we want to exclude anomalies from our data, but there are many situations where we don’t care about the “normal data” and are mostly interested in characterizing anomalies.   Supervised and Unsupervised Anomaly Detection   If we have historical data with a labeled target indicating the presence of anomalies in the data, we can use any supervised learning algorithm to perform anomaly detection. This is no different from traditional predictive modeling with a labeled target, although since anomalies are generally rare in the data, we often need to perform event-based sampling to ensure that the predictive model captures the relationship between the inputs and the target.   Consider a manufacturing line example where we have historical data about conveyor belt failures where the manufacturing line gets clogged, and our product starts to break and fail. These failures are rare (let’s say 99.9% of the data corresponds to normal conveyor belt operation, and 0.1% of the data corresponds to the failure mode), but we plan to increase profitability by predicting these conveyor belt failures before the product starts to break and fail. We could try to fit a predictive model like a neural network or a random forest to predict failure on the historical data, but our model would be 99.9% accurate in the training sample by predicting that failure never occurs. This is a generically good accuracy score for a predictive model, but in this case the model is useless because it fails to identify a single anomaly. A better approach would be to take an event-based sample of the data, so we select all 0.1% of failure cases and then randomly sample the remaining 99.9% of non-failure cases to get another 0.9% of the data. Overall, we only use 1% of the original data, but now anomalies are 10% of the sample and there is a much better chance that the model can capture the relationship between the inputs and the anomalies. The key idea is that our goals is to identify anomalies, so all the extra data about the non-failure cases does not help much in achieving our goal.   Predictive modeling algorithms are powerful tools for capturing relationships between inputs and targets, but when performing anomaly detection, we don’t always have labeled data and it isn’t necessarily feasible to collect this kind of data. In the previous example we had historical data about manufacturing line failures, but if we didn’t have this kind of data we probably don’t want to go out and break the manufacturing line to collect nice examples of anomalies. Without a labeled target, we must use unsupervised approaches to anomaly detection. The fundamental idea with unsupervised anomaly detection is to build a model that characterizes the normal data (through clustering, model building approaches, or semi-manual approaches), and then classify anything that deviates too much from the normal data as an anomaly. This approach is especially convenient for two related reasons: We can build anomaly detection models using normal operational data collected from functional systems or devices without the need to identify and label target information. Anomalies can be detected even if they are brand new events that have never been seen before. We can find anomalies that we were not necessarily “expecting” based on past data. It’s difficult to characterize all the ways that a system can break or all the methods that malign actors use to perpetrate fraud and malfeasance. On the other hand, it is usually easy to collect data and characterize the normal operation of a system. Unsupervised anomaly detection methods model the normal operational state and allow us to identify unusual behavior as anomalous even when this behavior is brand new.   Training and Deploying Anomaly Detection Models and Evaluating Performance   For supervised anomaly detection models, we must evaluate model performance on a validation sample that was not used to train the model. This means that we must first partition our historical data (which in this case consists of inputs and a target variable) into a training sample for building the model and a validation sample for evaluating model performance. Once we have built the model using the training sample, we score the validation sample to calculate model performance metrics. Standard assessment measures for binary targets like misclassification, accuracy, and area under the ROC curve can be useful for assessing model performance, but it is important to remember our modeling goal, which is to detect anomalies. As mentioned earlier, high accuracy is not useful in the model if we are not identifying anomalies with the model. Two assessment measures that are particularly useful for anomaly detection are precision and recall:   Precision = true positives / (true positives + false positives) Precision is a measure of how well our model distinguishes between anomalies and non-anomalies. It is the percentage of our “detected anomalies” (according to the model) that were labeled as anomalies in the data. Recall = true positives / (true positives + false negatives) Recall is a measure of our model’s ability to correctly identify anomalies. It is the percentage of total anomalies in the data that were identified by our model. High values of precision mean that most of the detected anomalies were truly anomalous, while high values of recall mean that our model doesn’t miss many anomalies in the data. Ideally, we would want high values of both precision and recall, but in practice we often must make a trade-off between models that detect a lot of anomalies and models that do a good job avoiding false positives. Deciding how to make this trade-off usually depends on the problem context, especially the cost associated with missing anomalies (too many false negatives) versus identifying more anomalies than exist (too many false positives).   Let’s go back to our conveyor belt manufacturing line example, in this case we want to use anomaly detection to identify conveyor belt failures before they cause a clog and destroy our products. For each anomaly we detect we will temporarily shut down the conveyor belt and manually inspect the products to ensure that they won’t get clogged. If we fail to detect an anomaly, the products will get clogged and start to break, requiring manual intervention to fix the clog and remove the waste. In this case false positives would waste money by forcing us to turn off the conveyor belt and manually inspect it for no reason, while false negatives would lead to clogged products and waste. The relative cost of these two activities (manual inspection and downtime vs. broken product and waste) will determine whether we care more about maximizing precision or recall. More generally the business costs associated with false positives and false negatives can be used to help determine how sensitive the model should be when detecting anomalies. If we are detecting anomalies to avoid catastrophic failure, we might be willing to accept a precision of 0.5 (so 50% of the detected anomalies are just normal cases) to achieve a recall of 0.99 (we capture 99% of all anomalies with the model).   If we have a labeled target and are using supervised learning techniques, we can just calculate the precision and recall on the validation data to evaluate model performance. If we are using unsupervised anomaly detection techniques, we will just have normal operating data with no labeled anomalies. In this case we need to deploy the model on real data or find examples of anomalies to see how the model performs. We can still calculate precision and recall, but only if we have new data containing known anomalies that we can use to judge model performance.     The main advantage of the unsupervised approach is that we don’t need labeled data to train the model, enabling us to build anomaly detection models in situations where it is simply impossible to build predictive models. If we do have labeled data, we can still perform honest assessment on these unsupervised models by seeing how they work on validation data containing anomalies. Training these models will involve only the non-anomalous data, so if we have labeled data our validation sample will include all the anomalies, and to do a good job with honest assessment it should also include some normal data.   Demonstration: Comparing Supervised and Unsupervised Methods for Anomaly Detection   Let’s make this discussion more tangible by exploring a couple of examples, one using a gradient boosted decision tree model (a popular machine learning algorithm) to perform supervised anomaly detection, and one using a support vector data description model (a popular unsupervised anomaly detection algorithm). We will skip the details of the algorithms for now and focus on the data setup for training and evaluating model performance. For this example, we will use the ionosphere data from the UCI machine learning repository, it contains antennae readings from radar signals that were used to probe the ionosphere to map its structure. Some of the radar signals hit the ionosphere and return to the antennae, providing information about the structure of the ionosphere, but some of the radar signals pass right through the ionosphere, revealing nothing. These useless radar signals are the anomalies we want to detect. This dataset does have labeled target called class, indicating whether the radar signals was anomalous. The inputs for this dataset are 32 variables, var000 through var0031, corresponding to the real and imaginary components of 16 complex electromagnetic waves (radar signals) measured by 16 different antennae. We try to use these radar variables to classify signals as anomalous or normal. cas; caslib _all_ assign; proc casutil; load file="/greenmonthly-export/ssemonthly/homes/Ari.Zitin@sas.com/ionosphere.csv" casout="ionosphere" replace; quit; proc freq data=casuser.ionosphere; table class; run; This is a bit of a toy example since this dataset only has 351 observations, 126 of which are anomalies. So already we have a good balance of normal observations and anomalies, and of course we have a labeled target variable indicating anomalous cases. proc partition data=casuser.ionosphere samppct=20 partind; by class; output out=casuser.ionosphere; run; We start by partitioning the data into a training sample and a validation sample. This is possible because we have a labeled target variable (class), and it is necessary to perform honest assessment when building the supervised predictive model. proc gradboost data=casuser.ionosphere(where=(_PartInd_ = 0)) outmodel=casuser.ionosphere_gradboost_model seed=919; input var_0000-var_0031 / level=interval; target class / level=nominal; run; We fit the gradient boosting model using the GRADBOOST procedure, using only the training data to fit the model. There are some hyperparameters we can choose for gradient boosting models, and we can even autotune these models. For this example, we will stick with the default settings on both the gradient boosting algorithm and the support vector data description algorithm, but when building models to deploy it is important to tune the hyperparameters to improve model performance. We output a model table that we use for scoring the validation data separately. proc gradboost data=casuser.ionosphere(where=(_PartInd_ = 1)) inmodel=casuser.ionosphere_gradboost_model; output out=casuser.gb_valid_scored copyvar=class; run; This second run of the GRADBOOST procedure scores the validation data using the trained model and then outputs the scored data into an in-memory CAS table. We will come back to this table and calculate precision and recall after we fit and score the unsupervised model. proc svdd data=casuser.ionosphere(where=(_PartInd_ = 0 and class=0)); id id; input var_0000-var_0031 / level=interval; kernel rbf / bw=mean; savestate rstore=casuser.ionosphere_svdd_astore; run; We use the SVDD procedure to fit the unsupervised support vector data description model. At the most basic level this model tries to draw a boundary tightly around the normal data (in this example in the 32-dimensional input space), and then all observations inside the boundary are normal and all observation outside the boundary are anomalous. Notice that when training this model, we use the training data (_PartInd_ = 0), but we also only trained on the normal data (class = 0). This is an situation where we could train the model even without a labeled target, although we would not be able to evaluate performance without at least some examples of the anomalies. Here we do have labeled target data, but we will only use it for model performance evaluation. We output a model ASTORE file that we can use for scoring the validation data separately. proc astore; score data=casuser.ionosphere(where=(_PartInd_ = 1)) rstore=casuser.ionosphere_svdd_astore out=casuser.svdd_valid_scored copyvar=class; run; We use the ASTORE procedure to score the validation data using the unsupervised SVDD model. In this case we include both the normal data and the anomaly data to see how the model does at distinguishing the anomalies from the normal data. This contrasts with the training process, where we only provided the normal data and did not include any examples of anomalies. Now that we have scored data for both the supervised and unsupervised modeling approaches, we can calculate precision and recall for each model. ods output ROCInfo=gb_roc; proc assess data=casuser.gb_valid_scored; var P_class1; target class / event="1" level=nominal; run; data work.gb_stats; set work.gb_roc; where (CutOff < 0.51 and Cutoff > 0.50); precision = TP / (TP+FP); recall = TP / (TP+FN); keep precision recall; run; data casuser.svdd_valid_scored; set casuser.svdd_valid_scored; if _SVDDSCORE_ = -1 then PredClass = 0; if _SVDDSCORE_ = 1 then PredClass = 1; run; ods output ROCInfo=svdd_roc; proc assess data=casuser.svdd_valid_scored; var PredClass; target class / event="1" level=nominal; run; data work.svdd_stats; set work.svdd_roc; where (CutOff < 0.51 and Cutoff > 0.50); precision = TP / (TP+FP); recall = TP / (TP+FN); keep precision recall; run; title 'Supervised Gradient Boosting Anomaly Detection'; proc print data=work.gb_stats noobs; run; title 'Unsupervised Support Vector Data Description Anomaly Detection'; proc print data=work.svdd_stats noobs; run; For each modeling approach we use the ASSESS procedure to calculate the confusion matrix across a range of probability cutoffs. For the SVDD model we don’t have any predicted probabilities (just predictions about whether we have anomalies), but we can still use PROC ASSESS to calculate the confusion matrix. Once we have the confusion matrix for the 0.5 probability cutoff we can calculate precision and recall for each model (the value of this cutoff matters and changing it could affect performance for the gradient boosting model but doesn’t matter for the SVDD model). We get the following results for the two models: Gradient Boosting Anomaly Detection Model: Precision = 0.96 Recall = 0.92 Support Vector Data Description Anomaly Detection Model: Precision = 0.49 Recall = 1.00   The gradient boosting model seems to outperform the SVDD model overall, although the SVDD model does capture all the anomalies in the data. The gradient boosting model does a better job of distinguishing between the anomalies and the normal data while the SVDD model is more sensitive and detects more anomalies. We could change the probability cutoff for the gradient boosting model to make it more sensitive, increasing recall at the cost of precision. The supervised learning approach ends up yielding a better model overall, but we can’t always choose this approach since we don’t always have labeled examples of anomalies. The advantage of the SVDD model is that we can still make progress in detecting anomalies even when we haven’t collected data about those anomalies. If we have a new conveyor belt on a manufacturing line and it has never broken, we don’t want to have to break it in a variety of ways to collect data on various kinds of anomalous failures. Instead, we can use an unsupervised anomaly detection method to look for deviations from the normal behavior of the new conveyor belt.   To build your own supervised and unsupervised anomaly detection models using SAS Viya check out the documentation below on the procedures used in this blog.   References/Useful Links: Ionosphere data set used in this example PROC GRADBOOST in SAS Viya PROC SVDD in SAS Viya ... View more

Watch this Ask the Expert session to learn how to build machine learning pipelines in SAS Viya.    Watch the Webinar   You will learn: How to use SAS Model Studio to build machine learning pipelines to improve your predictive modeling workflow. How nonlinear models, such as tree-based models, can improve model accuracy as compared to traditional linear models. How to add interpretability plots to black box models and how to use these plots to understand and improve predictive models.   The questions from the Q&A segment held at the end of the webinar are listed below and the slides from the webinar are attached.   Q&A Can you use a regression tree for a continuous dependent variable? Regression trees are suitable for a continuous dependent variable. In fact, if you build a decision tree model on a continuous target (the machine learning lingo for a dependent variable) you will be using a regression tree.   Does this have an open source node as does EM and can you utilize PROC Python within it? We do have the open source code node, and you don't even need to use PROC PYTHON. You can just connect to a Python install and execute Python (or R) code from a Model Studio node. This node can be used for data preprocessing or for machine learning, but does require data transfer (from the distributed computing environment to the Python environment)   Do you need to bring the data partitioning node to the pipeline? If it is done automatically, could you please show where it is done? The partition is done when the project is created in Model Studio. You have the option to turn it off, but you don't have to include a node for it (like in EM). This can be found in Model Studio under the Project Settings -> Partition Settings option.   Did your dataset of ~50,000 records include both training and validation data? Or did you upload two (or more) datasets? I uploaded one dataset and it was split automatically by the graphical interface. There is an option to manage the partition settings under Advanced -> Partition Settings while you are naming and creating the project. You can also access these settings from inside an already created project under Project Settings.   In the open source nodes, can you use native code or do you need to run R under Proc IML and Python under PROC Python? You can execute native Python or R code!   At the end of the class, can we go back to view a few charts to show us what may be signs of overfitting? (So we know what to watch out for.) If you see training performance is good but validation is bad, that is one of the best ways to tell.   My company is entertaining the idea of upgrading to Viya - will all the features be available, or do we have to sign (pay) for some of the features? We have separate products like base, EG, EM, etc. I'm so excited to hear you are thinking about moving to SAS Viya. It would be best to follow up with your SAS Account Executive and they can walk you through how the SAS Viya products translate to the SAS 9 products.   Do you have a Model Interpretability option in SAS EM? (I am still an EM user) There are macros. where you can build partial dependence semi manually. Now the macros will do it a little bit for you in SAS Enterprise Miner, but we don't have it included in the graphical interface in an easy way like we do in Model Studio. That is one of the advantages of upgrading. Here is a link to a paper with code for this in EM: https://www.sas.com/content/dam/SAS/support/en/sas-global-forum-proceedings/2018/1950-2018.pdf   Can these models automatically look for interactions in features, for example, handset age and customer age? The logistic regression model will not automatically look for interactions because it just creates the linear model. But the nonlinear models are fundamentally looking at interactions with how they're fitting the data. Basically, they look for nonlinear relationships between the inputs and the target, but they're not going to create an interaction variable and report out information about it. It's going to be mixed up behind the scenes in the complex nonlinear function that relates the inputs to the target. With decision trees, the nonlinear function is all those splits – those if then sequences. But for something like a neural network, it's going to be some complex nonlinear function of all the input variables.   What will be considered a threshold for ROC or is it always based on the training ROC value? A threshold like the worst RC, which is 0.5 for the area under the ROC curve is a random guessing model, so you want to be above that. In my mind, a value of around 0.6 is a weak model, whereas a value above 0.7 is a strong model. I would compare it to the training data to make sure we're not overfitting, but ultimately you want a model that predicts. Well, it does a good job making predictions on both the training and the validation data, so I usually want to make sure my ROC is above about 0.7. But for absolute scale, I generally like to look at the misclassification rate because that's easy to interpret in the business context. We want our misclassification rate to be less than 10%, for example, so that's kind of how I would set a standard rather than comparing it to the training data.   For local interpretability plots, is there a way to select which observation you are looking at? Yes, you can specify up to 5 observations. You can also specify a key variable in the Data tab and then use that variable to select specific observations from the data.   When the target is imbalanced where one outcome is much less likely to occur than it's complement (ex. cancer diagnosis), is there an option that improves prediction in forests? I know the randomForestSRC package includes the imbalanced function that samples the rarer event more frequently. Does SAS have a similar option available? Generally, you don’t need to balance classes unless the outcome is very rare (<10%), but there is the option to do event-based sampling in Model Studio. You would modify these settings in the project settings at the same time and place you choose the partition settings (before anything is ran in the pipeline).   Do you have guidelines for sample size for Trees? In general, you want to have enough training data to build a sensible model. This really depends on the kind of data and the modeling goal, but a good rule of thumb is to have at least 1000 observations in the training data for the tree-based models, although more observations will lead to a better model. A smaller sample size for the training data means a higher chance of overfitting, so that is the main thing to watch out for with small sample sizes. Also, keep in mind that the number of columns will play a large role in how many observations you need – as ideally you need multiple representations of every potential combination of columns.   How do we use custom function for auto-tuning? What if I wanted to select the hyperparameters based on Train and Test Gini? The autotuning algorithm allows you to select the assessment metric for the autotuning, which means you can choose how you want it to select hyperparameters. This can be found under Autotuning -> General Options in the node options. You can choose how the algorithm splits the training data and the objective function for autotuning.   How are the cut-offs in the nodes for the splits determined? Interval variables are binned so that they are not continuous. In the demo, Ari showed that the “Number of bins” was set to 100. This means that every numeric/interval variable is binned into 100 bins (can be set to quantile binning or bucket binning). Then, when the tree is determining splits it runs through all of the variables (binned numeric/interval and categorical) to determine what the best split would be and chooses that.   Recommended Resources SAS Tutorial: Machine Learning Fundamentals SAS Tutorial: Machine Learning: A Coding Example in SAS SAS Tutorial: How to Choose a Machine Learning Algorithm SAS Tutorial: Interpreting Machine Learning Models in SAS Moving from SAS®9 to SAS® Viya® Please see additional resources in the attached slide deck (PDF).   Want more tips? Be sure to subscribe to the Ask the Expert board to receive follow up Q&A, slides and recordings from other SAS Ask the Expert webinars.   ... View more