Method for real time network traffic classification

[0001] 1. Field

[0002] The present invention relates to computers and computer networks. More particularly, the present invention relates to network traffic classification method.

[0003] 2. Description of Related Art

[0004] Over the past few years, the Internet infrastructure has become a critical part of the global communications fabric. Emergence of new applications and protocols, such as voice over Internet Protocol (VoIP), pear-to-pear (P2P), and video conference increase the complexity of Internet. Accompanied with this trend is the increasing demand on more reliable and secure service. As a result, quality of service (QoS) becomes a more and more important issue. For example, transmission of real-time voice and video traffic has high bandwidth, low delay, and low packet loss rate requirements. Current best-effort Internet does not provide direct QoS guarantee to these real-time applications. At the same time, enterprise network administrators may want to restrict network bandwidth used by unauthorized VoIP or P2P applications. Hence, ability to guarantee or limit the QoS of specific types of network traffic may be desired in various different scenarios. Differentiated Services (DiffServ) is a method designed to guarantee different levels of QoS for different types of network traffic. It is achieved by configuring a “type of service” (TOS) field in a network data packet. For example, the DiffServ code point (DSCP) field in the IP header may be set according to the type of the network traffic for it to receive appropriate level of priority. Unfortunately, such design depends highly on network protocols, especially proprietary protocols, for observing DiffServ regulations. For example, in the worst case, the DiffServe method may be defeated if multiple protocols all set DSCP of their corresponding data packets to the highest number. Therefore, there is a need to classify network traffic flows in real time instead of relying on any tags in packet headers.

[0005] In general, in one aspect, the present invention involves classifying a network traffic flow by providing a first subspace, corresponding to a first classification, from a first training traffic flow; representing the network traffic in a stochastic process model; extracting a first power spectral density (PSD) feature vector from the network traffic flow according to a spectral analysis of the stochastic process model; measuring a first similarity of the first PSD feature vector with the first subspace according to a first similarity metric; and identifying a first classification of the network traffic flow according to the first similarity. Preferably, the present invention also includes providing a second subspace, corresponding to a second classification, from a second training traffic flow; measuring a second similarity of the first PSD feature vector with the second subspace according to a second similarity metric; and identifying a second classification of the network traffic flow according to the second similarity. Preferably, providing the first subspace includes representing the first training traffic in another stochastic process model; extracting a plurality of PSD feature vectors from the first training traffic flow according to the spectral analysis of the other stochastic process model; composing the first subspace using the plurality of PSD feature vectors; decomposing the first subspace into one or more segments; and identifying one or more bases each representing a structure of a segment. Preferably, the first subspace is decomposed according to minimizing a coding length of the one or more segments, the first similarity is measured by a normalized distance from the first PSD feature vector to the first subspace, and the normalized distance is defined according to the one or more bases of the first subspace. Preferably, the stochastic process model comprises a plurality of parameter set each corresponding to a packet of the network traffic flow, and each parameter set comprises a packet size and a packet arrival time associated with the corresponding packet. Preferably, the stochastic process model further includes a continuous time stochastic process; extracting the first power spectral density (PSD) feature vector includes processing the continuous time stochastic process through a low pass filter; sampling the low pass filtered continuous time stochastic process at a sampling frequency to produce a discrete time sequence; and extracting the PSD feature vector from a plurality of PSDs estimated from the discrete time sequence. Preferably, the sampling frequency corresponds to a period of substantially 0.5 milliseconds. Optionally, the first or second classification is one from the list consisting of audio traffic and video traffic. Optionally, the present invention also includes extracting a first sub flow of the network traffic flow according to the packet size; comparing a variance of the packet size to a predetermined threshold; and identifying a third classification of the network traffic according to the comparison. Optionally, the third classification is file traffic. Other aspects and advantages of the invention will be apparent from the following description and the appended claims.

[0017] A network traffic flow may be classified into multiple classifications, e.g., types or categories such as AUDIO (e.g., VoIP traffic flow), VIDEO (e.g., video conference traffic flow), FILE (e.g., file transfer traffic flow), UNKNOWN (e.g., unknown traffic flow), or other categories. Here a flow may be uniquely determined by 5 tuples in Transmission Control Protocol (TCP) and User Datagram Protocol (UDP), i.e., source IP address, destination IP address, source port, destination port, and protocol number. However, this classification scheme may be inappropriate because of the trend for session layer to reuse transport layer resources. For example, Skype, a proprietary peer-to-peer VoIP network, may be used by two users to communicate with each other while they may talk, type instant messages (IM), and transfer files at the same time. Skype may transmit different types of traffic over the same transport layer connection. Other applications such as MSN messenger (MSN), Google Talk (GTalk), or other instant messaging application may also allow session layer to reuse transport layer resources. In these examples, a network traffic flow may carry multiple types of traffic and may be referred to as “hybrid flow”. For example, a hybrid flow may be represented by a hypothesis vector {right arrow over ()}=<₀, ₁, ₂>, where

[0018] $\begin{array}{cc}{0}_{}=\{\begin{array}{cc}0& \mathrm{if}\mathrm{the}\mathrm{flow}\mathrm{carries}\mathrm{no}\mathrm{FILE}\mathrm{traffic}\\ 1& \mathrm{if}\mathrm{the}\mathrm{flow}\mathrm{carries}\mathrm{FILE}\mathrm{traffic}\end{array},& \left(1\right)\\ 1_{}=\{\begin{array}{cc}0& \mathrm{if}\mathrm{the}\mathrm{flow}\mathrm{carries}\mathrm{no}\mathrm{AUDIO}\mathrm{traffic}\\ 1& \mathrm{if}\mathrm{the}\mathrm{flow}\mathrm{carries}\mathrm{AUDIO}\mathrm{traffic}\end{array},& \left(2\right)\\ 2_{}=\{\begin{array}{cc}0& \mathrm{if}\mathrm{the}\mathrm{flow}\mathrm{carries}\mathrm{no}\mathrm{VIDEO}\mathrm{traffic}\\ 1& \mathrm{if}\mathrm{the}\mathrm{flow}\mathrm{carries}\mathrm{VIDEO}\mathrm{traffic}\end{array}.& \left(3\right)\end{array}$
In some examples, the hypothesis vector {right arrow over ()}=<₀, ₁, ₂> may represent a true state of a hybrid flow and {right arrow over ()}=<{circumflex over ()}₀, {circumflex over ()}₁, {circumflex over ()}₂> may represent an estimated hypothesis vector of the hybrid flow. Although only three parameters ₀, ₁, and ₂ are shown in the three dimensional hypothesis vector here, different number of parameters may be necessary and the hypothesis vector may have higher or lower dimension than three in other examples. In another example, different types of network traffic flow may be included in the hybrid flow and the hypothesis vector may represent network traffic flows different from the AUDIO, VIDIO, or FILE.

[0019] A network traffic flow may be modeled by as a stochastic process. For example, the network traffic flow may be defined as a set of data packets, each data packet being represented by data packet features, such as packet size and packet arrival time. Thus, a network traffic flow may be represented by
={_i, _i; i=1, . . . , I},  (4)
where _i and _i denote the packet size and relative packet arrival time, respectively, of the ith packet in the flow with total number of I packets. In one example, only relative arrival time may be considered in which case ₁ arrival is always 0. In other examples, there may be different number of data packet features and the data packet features may be different than described above.

[0020] FIG. 1A shows a system architecture diagram in accordance with one embodiment of the present invention. As shown in FIG. 1, the system architecture 100 includes classification module 122 and the training module 121. The classification module 122 includes network traffic flow 101, sub flow extraction module 102, PSD feature extraction modules 103, AUDIO/VIDEO classifier module 104, FILE classifier module 105, ₀ parameter 106, ₁ parameter 107, and ₂ parameter 108. The training module 121 includes AUDIO training traffic flow 109, VIDEO training traffic flow 112, PSD feature extraction modules 110 and 113, subspace identification modules 111 and 114, AUDIO subspace 115, and VIDEO subspace 116. There are many ways to implement sub flow extraction, for example, in FIG. 1A the network traffic flow 101 is split into two sub-flows _L and _S by the sub flow extraction module 102. The sub flow _L is composed of packets larger than a threshold θ_P and _S is composed of smaller packets. In one example, the FILE classifier module 105 determines the ₀ parameter 106 based on the sub flow _S. In another example, the sub flow _L is processed by the PSD feature extraction modules 103 and AUDIO/VIDEO classifier module 104 to determine the ₁ parameter 107 and the ₂ parameter 108. There are many ways to implement AUDIO/VIDEO classification, for example, in FIG. 1A the AUDIO/VIDEO classifier module 104 classifies the sub flow _L according to the AUDIO subspace 115 and VIDEO subspace 116. In an example shown in FIG. 1A, the AUDIO subspace 115 is generated from the PSD feature extraction module 110 and subspace identification module 111 according to the AUDIO training traffic flow 109 while the VIDIO subspace 116 is generated from the PSD feature extraction module 113 and subspace identification module 114 according to the VIDIO training traffic flow 112. Although only two sub flows _L and _S are shown in the example in FIG. 1A here, there may be more or less number of sub flows in other examples.

[0021] FIG. 2 shows a flow chart of a method for real time classification of network traffic flow in accordance with one embodiment of the present invention. Initially, a sub flow of a network traffic flow is extracted according to packet size of data packets in the network traffic flow (step 210). For example, file transfer traffic flow has two common properties. First, although it is not specified in network protocols, almost all file transfer applications are implemented in such a way that they partition a file into equal-sized segments and transfer one segment at a time. Therefore, a FILE traffic flow may be composed of equal-sized packets, except for the last packet, whose size is usually smaller because file size may not be a multiple of segment size. Secondly, packet size of FILE traffic flow is larger than that of AUDIO or VIDEO traffic flow. The larger the file segment is, the fewer the number of packets may be required to transmit a file, which translates to less overhead contributed by packet header. To avoid fragmentation, software developers may also want to choose packet size smaller than maximum transmission unit (MTU). The most common MTU value of Ethernet, and hence the majority of the Internet, is 1500 bytes. As a result, packet size of FILE traffic flow may commonly be chosen to be close to 1500 bytes. For example, Skype uses 1397 bytes and MSN uses 1448 bytes for each packet in a FILE traffic flow. On the other hand, packet size in AUDIO or VIDEO traffic flow may be smaller. By examining Skype, MSN, and GTalk, it may be observed that VoIP packet size is nominally from 50 bytes to 200 bytes, and video conference packet size is nominally from 50 bytes to 600 bytes. Therefore, based on the description above, packet size in a network traffic flow may be used as criteria to separate FILE traffic flow from AUDIO or VIDEO traffic flow within a hybrid flow. Continue to the discussion of FIG. 2, a variance of packet size in a network traffic flow is compared to a predetermined threshold in step 203. For example, once a sub flow _L is extracted in step 201, the variance of packet size in the sub flow _L may be estimated, for example by the FILE classifier module 105 shown in FIG. 1A. As mentioned above, FILE traffic flow is commonly composed of equal-sized packets, except for the last packet. Thus, variance of packet size in a FILE traffic flow may be close to zero. Therefore, a network traffic flow may be classified as FILE traffic flow according to the comparison in step 203 (step 205). In one example, in step 203 a threshold may be set as θ₉₄ and in step 205 the FILE classifier module 105 may set {circumflex over ()}₀ such that

[0022] $\begin{array}{cc}{\stackrel{}{^}}_0=\{\begin{array}{cc}0& \frac{1}{I-1}\sum _{i=1}^I{\left(i_{}-\stackrel{}{\_}\right)}^2>{\theta }_{\sigma }\\ 1& \frac{1}{I-1}\sum _{i=1}^I{\left(i_{}-\stackrel{}{\_}\right)}^2\le {\theta }_{\sigma }\end{array},\mathrm{where}\stackrel{}{\_}=\frac{1}{I}\sum _{i=1}^I.& \left(5\right)\end{array}$

[0023] FIG. 3 shows another flow chart of a method for real time classification of network traffic flow in accordance with one embodiment of the present invention. Comparing to the example shown in reference to FIG. 2 above, classification of AUDIO or VIDEO traffic flows is more complicated. Generally speaking, to solve a pattern classification problem, one defines a feature vector, which is a summary (also known as a statistic) of raw bit stream in a network traffic flow. During a training phase, the method may extract training feature vectors from training traffic flows for each category and maintains them or their statistics, such as a subspace as described in reference to FIG. 4A below (step 301). During a classification phase, a network traffic flow, or a sub flow extracted from a network traffic flow, may be first represented as a stochastic process as described above (step 303). The method may extract PSD feature vectors using the stochastic process model (step 305). A general classifier, e.g., the AUDIO/VIDIO classifier module 104 shown in FIG. 1A, then makes decision by measuring similarity between the feature vector of the network traffic flow being classified and each of the multiple training feature vectors (step 307). The classification is then determined according to predefined criterion (step 309).

[0024] Referring to step 305, there are many ways to define the mapping from raw bit stream in a network traffic flow to feature vector. In one example, power spectral density (PSD) feature vector, which exploits regularities in AUDIO and VIDEO traffic, may be extracted by PSD feature extraction modules 110, 113, and 103 shown in FIG. 1A. More details are described in reference to FIGS. 1C and 4 below.

[0025] Referring to step 307, there are also many ways to define the similarity metric. Many similarity metrics have been used in the prior art, for example, cost function, Euclidean distance of nearest-neighbor (1-NN) or K-nearest-neighbor (K-NN), and others. Generally speaking, no similarity metric is guaranteed to be the best for all applications. For example, certain cost function is applicable only when likelihood probabilities are well estimated, which requires the number of training traffic flows to be much larger than the number of dimensions in a feature vector. Therefore, it is not suitable for classification based on high dimensional feature vector, such as the PSD feature vector as described below in reference to FIG. 1B. Euclidean distance of 1-NN and K-NN are optimal only when the assumption is true that data, i.e., sampled data corresponding to network traffic flow, of same category are clustered together. This may not always be true. In an example, by studying the data set, a set of sampled data corresponding to a set of network traffic flows, a similarity metric may be based on the normalized distance from feature vector of the network traffic flow being classified to a subspace obtained during training phase, as described in step 301 above. The idea is the following. Suppose the feature vectors belong to a space ^M of dimensions M. During the training phase, low (i.e., less than M) dimensional subspace structures of training data may be identified from the high dimensional space ^M. Here, the low dimensional subspace is a portion, or subset, of the high dimensional space ^M. This may be identified by the subspace identification modules 111 and 114 as shown in FIG. 1A. More details are described in reference to FIG. 5 below.

[0026] Referring back to step 307, during the classification phase, the normalized distances may be calculated between the feature vector of the network traffic flow being classified and the identified subspaces. In one example, in each category, such as AUDIO traffic flow and VIDEO traffic flow, the normalized distance may be defined as the minimum distance. For example, by comparing the normalized distances, to the identified subspace associated with the AUDIO traffic flow category and the VIDEO traffic flow category respectively, with their corresponding predefined thresholds, the network traffic flow may be classified as AUDIO traffic flow, VIDEO traffic flow, or UNKNOWN traffic flow. This may be performed by the AUDIO/VIDEO classifier module 104 as shown in FIG. 1A.

[0027] AUDIO and VIDEO protocols may specify a constant interval between two consecutive packets at the transmission side. This constant interval may be referred to as inter-packet delay (IPD). For example, table I lists several speech codec standards and the corresponding specified.

[0028]

TABLE I
Commonly used speech codec and their specifications
	Standard	Codec Method	Inter-Packet Delay (ms)
	G.711 [28]	PCM	.125
	G.726 [29]	ADPCM	.125
	G.728 [30]	LD-CELP	.625
	G.729 [31]	CS-ACELP	10
	G.729A [31]	CS-ACELP	10
	G.723.1 [32]	MP-MLQ	30
	G.723.1 [32]	ACELP	30

[0029] Due to random propagation delay of each packet, the inter-arrival time (IAT) between consecutive packets arriving at a network device, e.g., a router where network traffic features are measured, may also be random. This is known as jitter. However, the IATs may not vary too much if IPDs remains constant. This may constitute a regularity pattern residing in packet arrival time parameters of an AUDIO or VIDEO traffic flow. In addition, packet sizes of AUDIO or VIDEO traffic flows vary periodically. For example, most video coding schemes use two types of frames, i.e., intra frames (I-frame) and predicted frames (P-frame). An I-frame is a frame coded without reference to any other frame except itself. It serves as the starting point for a decoder to reconstruct the video stream. A P-frame may contain both image data and motion vector displacements and/or combinations of the two. Its decoding needs reference to previously decoded frames. Packets containing I-frames are in general larger than those containing P-frames. Usually, the number of P-frames between two consecutive I-frames is constant. Therefore, periodic variation of packet size may be observed. Voice streams may have similar phenomenon if linear prediction coding (LPC), e.g., code excited linear prediction (CELP) voice coder, is employed. Therefore, multiple regularity patterns may reside in packet size parameters of a network traffic flow. Further, such regularity pattern may remain invariable in a hybrid flow. FIG. 1B shows an exemplary network traffic flow in accordance with one embodiment of the present invention. FIG. 1B depicts a scenario when a network traffic flow carries multiple sessions, such as a VoIP session and a file transfer session. At the transmission side, AUDIO packets 131 are transmitted periodically whereas FILE packets 132 are transmitted more randomly. When the flow is captured, the high regularity of AUDIO packet arrival events still maintains, although packets of two sessions are interleaved with each other. Accordingly, by identifying regularities in network traffic flows, especially hybrid flows, AUDIO or VIDEO traffic flows may be classified. In one example, time domain regularity detection may be employed but it is not a straight-forward task. In other examples, PSD features, which contain information in frequency domain may be employed.

[0030] FIG. 1C shows a schematic block diagram of the PSD feature extraction module 150 in accordance with one embodiment of the present invention. The PSD feature extraction modules 103, 110, and 113, as shown in FIG. 1A, are essentially the same as of the PSD feature extraction module 150. Here, the PSD feature extraction module 150 includes preprocessing module 151, low pass filter (LPF) 152, sampling module 153, and PSD estimate module 154. The preprocessing module 151 converts a network traffic flow _S into a stochastic process model (t). The low pass filter (LPF) 152 converts stochastic process model P(t) into a low pass filtered output P_h (t). The sampling module 153 converts P_h (t) into discrete samples P_d(i). The PSD estimate module 154 converts the discrete samples P_d(i) into {right arrow over (ψ)}ε^M.

[0031] FIG. 4 shows yet another flow chart of a method for real time classification of network traffic flow in accordance with one embodiment of the present invention. Initially, a network traffic flow _S, which may be modeled in a stochastic process as described in equation (4) above, is converted to a continuous-time stochastic process below:

[0032] $\begin{array}{cc}\left(t\right)=\sum _{<,>\in }\delta \left(t-\right),& \left(6\right)\end{array}$
where δ(·) denotes the delta function (step 402). Because digital computers are more suitable to deal with discrete-time sequences than continuous-time processes, it is desirable to transform (t) to discrete-time sequence. This may be obtained by digital sampling. One critical concern over digital sampling is the aliasing problem. Suppose the sampling interval is T_s· then the sampling frequency

[0033] $F_s=\frac{1}{T_s}.$
By Nyquist-Shannon sampling theorem, the sampling frequency F_s should be greater than twice the maximum frequency component in the original continuous-time signal. Otherwise, aliasing happens. That is, if the original signal has a frequency component F′ such that

[0034] $F\text{'}>\frac{F_s}{2}$
then, after sampling, this frequency component is mapped to frequency F″ in the sampled sequence such that
F′≡F″mod F′_s,  (7)

[0035] $\begin{array}{cc}F\mathrm{\text{'}\text{'}}<\frac{F_s}{2}.& \left(8\right)\end{array}$
Once aliasing happens, information in high frequency conflicts with information in low frequency. This is not desired. The continuous-time stochastic process, shown in equation (6), is a summation of delta functions, whose spectra span the entire frequency domain. Therefore, it is necessary to reshape the spectrum of (t) to avoid aliasing when it is sampled at an interval T_s·. The step 404 may be accomplished by a LPF module 152, having a transfer function h_LPF(t), whose output may be denoted as _h(t), as shown in FIG. 1C and equation (9) below.

[0036] $\begin{array}{cc}{h}_{}\left(t\right)=\left(t\right)*h_{LPF}\left(t\right)=\sum _{<,>\in s}h\left(t-\right).& \left(9\right)\end{array}$
Then _h(t) may be sampled at an interval T_s· to produce a discrete-time sequence, as shown in equation (10) below (step 406).

[0037] $\begin{array}{cc}{d}_{}\left(i\right)=h_{}\left({\mathrm{iT}}_s\right)=\sum _{,}^{}h\left({\mathrm{iT}}_s-\right),& \left(10\right)\end{array}$
where i=1, . . . ,

[0038] $I_d=\frac{{\max }_{}}{T_s}+1,{\max }_{}$
is the arrival time of the last data packet in the network traffic flow _S·.

[0039] The sampling interval T_s· may not be arbitrarily chosen. If T_s· is too large, then F_s is too small such that useful information in high frequency may be lost. If T_s· is too small, then the discrete-time sequence length, I_d·, is too large, which translates to high complexity in PSD estimation. As described above in reference to FIG. 1B, one type of regularity pattern comes from packet arrival events in real-time network traffic flow. In one example,

[0040] $T_s<\frac{I{\mathrm{AT}}_{\min }}{2},$
where IAT_min stands for the minimum IAT in the network traffic flow, may be chosen to substantially retain all useful information. However, when multiple types of network traffic flows are interleaved in one hybrid flow, IATs are much smaller than those of a non-hybrid flow. In such case,

[0041] $\frac{I{\mathrm{AT}}_{\min }}{2}$
may be smaller than necessary. In another example, by studying common software supporting VoIP and video conference, such as Skype, MSN, and GTalk, T_s· of 0.5 milliseconds may be chosen to substantially extract all useful information for the purpose of AUDIO and VIDEO traffic flow classification.

[0042] After sampling, the PSD feature vector may be extracted by way of estimation from the discrete-time sequence (step 408). This may be accomplished by identifying regularity patterns in (t), or by identifying periodic components in the discrete-time sequence {_d(i)}_i=1^Id. In one example, identification of regularity pattern may be performed in time domain, which may be complex. In another example, PSD feature vector of {_d(i)}_i=1^Id may be estimated in frequency domain. For a second-order stationary sequence {y_i; iε}, its PSD feature vector may be defined as

[0043] $\begin{array}{cc}\psi \left(\varpi ;y\right)=\sum _{k=-\infty }^{\infty }r\left(k;y\right)e^{-j\varpi k},\varpi \in \left[-\pi ,\pi \right),& \left(11\right)\end{array}$
where {r(k;y); kε} is the auto-covariance sequence (ACS) of {y_i; iε}, i.e.,
r(k;y)=ξ[y(i)y*(i−k)].  (12)
The parameter ω in equation (11) may be restricted to be within [−π, π) because ψ( ω; y)=ψ( ω+2π; y). PSD feature vector may represent the energy distribution of the sequence {y_i; iε} in frequency domain. Regularity patterns in {y_i; iε} translate to dominant periodic components in its ACS, which may further map to peaks in its PSD feature vectors. This is because PSD feature vector is simply the discrete-time Fourier transform (DTFT) of ACS. According to equations (11) and (12), obtaining true PSD feature vectors requires infinitely long time series. In some examples, PSD feature vectors may be estimated from a finite portion of the sequence. There are many PSD estimation techniques, such as nonparametric methods, e.g., Welch method and filter-bank methods, and parametric methods, e.g., autoregressive (AR) model, moving average (MA) model, and autoregressive moving average (ARMA) model. In one example, AR model may be used to model signals whose PSD have sharp peaks, i.e., high regularity in time domain. In another example, Levinson-Durbin Algorithm (LDA) may be used to estimate the AR model for obtaining the PSD. For example, the discrete-time sequence {_d(i)}_i=1^Id may be assumed to be second-order stationary. Then its PSD may be represented as
ψ( ω;_d), ωε[−π, π).  (13)
Recall that {_d(i)}_i=1^Id are obtained by sampling a continuous-time signal _h(t) at time interval T_s·, thus, the PSD may be further represented in terms of real frequency as

[0044] $\begin{array}{cc}{\psi }_f\left(f;d_{}\right)=\psi \left(\frac{2\pi f}{F_s};d_{}\right),& \left(14\right)\end{array}$
where

[0045] $F_s=\frac{1}{T_s}.$
Equation (14) shows the relationship between periodic components in continuous-time stochastic process and shape of PSD in frequency domain.
ψ_f(f;_d) is a continuous function in frequency domain. In an example, a series of frequencies may be selected:

[0046] $\begin{array}{cc}0\le f_1<f_2<\dots <f_M\le \frac{F_s}{2},& \left(15\right)\end{array}$
and the PSD feature vector may be defined as
{right arrow over (ψ)}=[ψ_f(f₁;_d), ψ_f(f₂;_d), ψ_f(f_M;_d)]^T.  (16)
{right arrow over (ψ)}ε^M.

[0047] As shown in FIG. 1A, PSD feature extraction modules 103, 110, and 113 are used in both training and classification phases. The time complexity of PSD feature extraction module may be composed of two parts, i.e., time complexities of preprocessing and PSD estimation. From Equation (10), time complexity of preprocessing part is O (|_S|). The time complexity of the LDA is O(I_dlog(I_d)). As |_S|<<I_d·. The total time complexity may be represented substantially as:
O(I_dlog (I_d)).  (17)

[0048] FIG. 5 shows still another flow chart of a method for real time classification of network traffic flow in accordance with one embodiment of the present invention. Initially, a training traffic flow is represented in a stochastic process model (step 501). Then PSD feature vectors are extracted from the training traffic flow according to spectral analysis of the stochastic process model (step 503). A subspace may then be composed using the PSD feature vectors (step 505). The subspace may be decomposed into segments (step 507) and then subspace bases each representing a structure of a segment may be identified (step 509).

[0049] As shown in FIG. 1A, AUDIO and VIDEO training traffic flows may be processed during the training phase. After processing these raw data through PSD feature extraction modules 110 and 113, two sets of feature vectors may be obtained as:
1) Ψ⁽¹⁾Δ{{right arrow over (ψ)}₁(1), {right arrow over (ψ)}₁(2), . . . , {right arrow over (ψ)}₁ (N₁)},  (18)
which is obtained through processing AUDIO training traffic flows, where N₁ is the number of AUDIO training traffic flows; and
2) Ψ⁽²⁾Δ{{right arrow over (ψ)}₂(1), {right arrow over (ψ)}₂(2), . . . , {right arrow over (ψ)}₂ (N₂)},  (19)
which is obtained through VIDEO training traffic flows, where N₂ is the number of VIDEO training traffic flows. To facilitate further discussion, let us also regard Ψ⁽ⁱ⁾ as a M×N_i matrix, for i=1, 2, where each column is a feature vector. In other word, Ψ⁽ⁱ⁾ε^M×Ni.

[0050] Referring to step 505, there are many techniques to identify the low dimensional subspaces embedded in ^M for both Ψ₁ and Ψ₂, such as principal components analysis (PCA) and Metric Multidimensional Scaling (MDS), which identify linear structure, and ISOMAP and locally linear embedding (LLE), which identify non-linear structure. All these methods assume that data are embedded in one single low-dimensional subspace. This assumption may not always be true. For example, as different software uses different voice coding, the PSD feature vector of AUDIO traffic flow may be a random vector generated from a combination of models than a single model. In such case, it is more likely that there may be multiple subspaces in which the feature vectors are embedded. Same may hold true for feature vectors of VIDEO traffic flows. In some examples, the method may include 1) clustering training feature vectors into multiple groups, which may be referred to as subspace decomposition; and 2) identifying subspace structure in each group, which may be referred to as subspace bases identification.

[0051] Referring to step 507, the purpose of subspace decomposition is to partition the data set Ψ, which includes N vectors each known as a base:
Ψ={{right arrow over (ψ)}(1), {right arrow over (ψ)}(2), . . . , {right arrow over (ψ)}(N)}  (20)
into non-overlapping K subsets such that
Ψ=Ψ₁∪Ψ₂∪ . . . ∪Ψ_K.  (21)
For example, subspaces may be decomposed according to a minimum coding length criteria. The idea is to view data segmentation problem from the perspective of data coding/compression. Suppose one wants to find a coding scheme , which maps data in Ψε^M×N to a bit sequence. As all elements are real numbers, infinite long bit sequence is needed to decode without error. Hence, one has to specify the tolerable decoding error, ε to obtain a mapping with finite coding length, i.e.,
∥{right arrow over (ψ)}_n−⁻¹(({right arrow over (ψ)}_n))∥²≦ε², for ∀n=1, . . . , N.  (22)
Then the coding length of the coding scheme is a function
: ^M×N→₊.  (23)
The optimal partition should minimize coding length of the segmented data, i.e.,

[0052] $\begin{array}{cc}\begin{array}{c}\min \\ \Pi \end{array}\stackrel{}{^}\left(\Psi ;\Pi \right)=\begin{array}{c}\min \\ \Pi \end{array}\left\{\sum _{k=1}^{{\rm K}}\left({\Psi }_k\right)+\sum _{k=1}^{{\rm K}}{\Psi }_k\left[-{\log }_2\left(\frac{{\Psi }_k}{{\rm K}}\right)\right]\right\},& \left(24\right)\end{array}$
where II is the partition scheme. The first term in equation (24) is the summation of coding length of each group, and the second one is the number of bits needed to encoding membership of each item of Ψ in the K groups. The optimal partition may be achieved in the following way. Let the segmentation scheme be represented by the membership matrix,
π_kΔdiag ([π_1k, π_2k, . . . , π_Nk])ε^N×N,  (25)
where π_nk denotes the probability that vector {right arrow over (ψ)}(n) belongs to subset k such that

[0053] $\begin{array}{cc}\sum _{k=1}^{{\rm K}}{\pi }_{\mathrm{nk}}=1,\mathrm{for}\forall n=1,\dots ,N& \left(26\right)\end{array}$
and diag(.) denotes converting a vector to a diagonal matrix. It is known in the art that the coding length may be upper bounded that

[0054] $\begin{array}{cc}\stackrel{}{^}\left(\Psi ;\Pi \right)\le \sum _{k=1}^{{\rm K}}\left[\frac{\mathrm{tr}\left({\Pi }_k\right)+M}{2}{\log }_2\det \left(I+\frac{M}{\mathrm{tr}\left({\Pi }_k\right){\epsilon }^2}\Psi {\Pi }_k{\Psi }^T\right)\right]+\sum _{k=1}^{{\rm K}}\left[\mathrm{tr}\left({\Pi }_k\right)\left(-{\log }_2\frac{\mathrm{tr}\left({\Pi }_k\right)}{N}\right)\right]\stackrel{\Delta }{=}\stackrel{}{^}\left(\Psi ;\Pi \right),& \left(27\right)\end{array}$
where tr (.) denotes the trace of a matrix, and det (.) denotes matrix determinant. Combining equations (24) and (27), a minimax criterion may be obtained as:

[0055] $\begin{array}{cc}\hat{\Pi }=\arg \underset{\Pi }{\min }\left[\underset{}{\max }\stackrel{}{^}\left(\Psi ;\Pi \right)\right]=\arg \underset{\Pi }{\min }\Psi ;\Pi ).& \left(28\right)\end{array}$
Then for ∀ψ_nεΨ, ψ_nεΨ_k after segmentation if and only if

[0056] $\begin{array}{cc}\hat{k}=\arg \underset{k}{\min }{\hat{\pi }}_{\mathrm{nk}}.& \left(29\right)\end{array}$
There is no closed form solution for equation (28). It is known in the art that a pair-wise steepest descent method may be used to solve it. Using the above method, a partition of AUDIO feature vector set Ψ⁽¹⁾
Ψ⁽¹⁾=Ψ₁⁽¹⁾∪ . . . ∪Ψ_K1⁽¹⁾,  (30)
and a partition of VIDEO feature vector set Ψ⁽²⁾
Ψ⁽²⁾=Ψ₁⁽²⁾∪ . . . ∪Ψ_K2⁽²⁾.  (31)
may be obtained.

[0057] Referring to step 509, an exemplary algorithm to identify subspace bases for each segmentation,
{Ψ_k⁽ⁱ⁾; k=1, . . . , K_i, i=1,2},  (32)
may be shown as the function IdentifyBases below:

[0058]

1)	function [{right arrow over (μ)}, Û, {circumflex over (Σ)}, Ū, Σ] = IdentifyBases(Ψ ε M×N, δ)
2)	μ → = 1 → | Ψ | ⁢ ∑ ψ → ⁢ ∈ Ψ ⁢ ψ →
3)	Ψ = [{right arrow over (ψ)}1 − {right arrow over (μ)}, {right arrow over (ψ)}2 − {right arrow over (μ)}, . . . , {right arrow over (ψ)}|Ψ| − {right arrow over (μ)}]
4)	Do eigenvalue decomposition on Ψ ΨT such that
	Ψ ΨT = UΣUT,	(33)
	where U [{right arrow over (μ)}1, . . . , {right arrow over (μ)}M], Σ diag ([σ12, . . . , σM2]), and
	σ12 ≧ σ22 ≧ . . . ≧ σM2.
5)	= arg ⁢ ⁢ min ⁢ ∑ m = 1 ⁢ ⁢ σ 2 ⁡ ( m ) ≥ δ ⁢ ∑ m = 1 M ⁢ ⁢ σ 2 ⁡ ( m )
6)	Û = [{right arrow over (u)}1, {right arrow over (u)}2, . . . , {right arrow over (u)} −1]
7)	Ū = [{right arrow over (u)} , {right arrow over (u)} +1, . . . , {right arrow over (u)}M]
8)	{circumflex over (Σ)} = diag ([σ12, . . . , σ −12])
9)	Σ = diag ([σ 2, . . . , σM2])
10)	end function

The basic idea is to identify uncorrelated bases and choose those bases with dominant energy. As shown in equation (33), argument Ψ is the feature vector set of one segmentation and δ is a user defined parameter which specifies the percentage of energy retained, e.g., 90% or 95%. It returns 5 variables. {right arrow over (μ)} is the sample mean of all feature vectors. It is the origin of the identified subspace. Columns of {circumflex over ()} are the bases with dominant energy (i.e., high concentration of occurrence in the frequency domain distribution), whose corresponding variances are {circumflex over (Σ)}. These bases determine the identified low dimensional subspace spanned by Ψ. Columns of compose the null space of the previous subspace, whose corresponding variances are Σ. The last two outputs are required to calculate distance of a test feature vector to the subspace, which will be described in more details below. Applying the function IdentifyBases on all segmentation
[{right arrow over (μ)}_k⁽ⁱ⁾,{circumflex over ()}_k⁽ⁱ⁾, {circumflex over (Σ)}_k⁽ⁱ⁾,_k⁽ⁱ⁾, Σ_k⁽ⁱ⁾]=IdentifyBases(Ψ_k⁽ⁱ⁾)  (34)
for ∀k=1, . . . , K_i, ∀i=1, 2, may be obtained. These are the exemplary outputs of subspace identification modules 111 and 114, or the results of training phase AUDIO subspace 115 or VIDEO subspace 116 as shown in FIG. 1A.

[0059] As described above, an exemplary approach to identify subspaces spanned by PSD feature vectors of AUDIO or VIDEO training traffic flows. Specifically, parameters
[{right arrow over (μ)}_k⁽ⁱ⁾, {circumflex over ()}_k⁽ⁱ⁾, {circumflex over (Σ)}_k⁽ⁱ⁾, _k⁽ⁱ⁾, Σ_k⁽ⁱ⁾]  (35)
for ∀k=1, . . . , K_i, ∀i=1, 2, may be obtained for real time classification. Referring back to the classification phase 122 in FIG. 1A, given a network traffic flow , a sub flow _S may be extracted, by extracting small packets in the sub flow extraction module 102, i.e., packets smaller than θ_P, and passes it through the PSD feature extraction module 103 to generate PSD feature vector {right arrow over (ψ)}. This is the input to the AUDIO/VIDEO classifier module 104. For example, the AUDIO/VIDEO classifier module 104 works in the following way in step 307 as described in reference to FIG. 3. It first calculates normalized distances between {right arrow over (ψ)} and all subspaces of both categories. Then it chooses minimum distance to each category. By comparing the two distance values to two thresholds, θ_A and θ_V, respectively, decision is made. An exemplary procedure of the AUDIO/VIDEO classifier may be shown below:

[0060]

	1)	function [ 1, 2] = AUDIOVIDEOClassify ({right arrow over (ψ)}, θA, θV)
	2)	For ∀i = 1, 2, ∀k = 1, . . . , Ki
		dk(i) = NormalizedDistance ({right arrow over (ψ)}, {right arrow over (μ)}k(i), Ūk(i), Σk(i))
	3)	For ⁢ ⁢ ∀ i = 1 , 2 , ⁢ d t = min k ⁢ d k ( t )
	4)	if d1 < θA and d2 > θV
	5)	1 = 1, 2 = 0, i.e., AUDIO flow.
	6)	else if d1 > θA and d2 < θV
	7)	1 = 2 = 1, i.e., VIDEO flow.
	8)	else
	9)	1 = 2 = 0, i.e., neither AUDIO nor VIDEO.
	10)	end if
	11)	end function
	12)	function d = NormalizedDistance({right arrow over (ψ)}, {right arrow over (μ)}, Ū, Σ)
	13)	d = ({right arrow over (ψ)} − μ)T Ū Σ−1ŪT ({right arrow over (ψ)} − {right arrow over (μ)})
	14)	end function

In some examples, it may not able to distinguish a pure VIDEO traffic flow from a VIDEO+AUDIO hybrid flow. First, the difference of a VIDEO traffic flow and a VIDEO+AUDIO hybrid flow may be very small, because the piggyback technique is widely used such that AUDIO and VIDEO data may be encapsulated in the same packet. Secondly, it may be of little practical value to do that as it is very rare for two users watching each other without talking. Accordingly, in the function AUDIOVIDEOClassify, line 7, when {circumflex over ()}₂=1, the algorithm always sets {circumflex over ()}₁=1. In other examples, the function AUDIOVIDEOClassify may be implemented with a functionality to distinguish a pure VIDEO traffic flow from a VIDEO+AUDIO hybrid flow. The parameters ₁ and ₂ may be determined independly.

[0061] From lines 2 and 13 above, the time complexity of AUDIO and VIDEO classifier is
O ((K₁+K₂)M²).  (36)
Combining equations (17) and (36), the overall complexity of the real time classification scheme may be represented as
O (I_d log(I_d))+O((K₁+K₂)M²).  (37)

[0062] Experiment results of applying the exemplary system architecture of FIG. 1A to real time network traffic classification may be summarized here. For example, two sets of experiments are conducted on network traffic flows generated by Skype with the results shown in FIGS. 6 and 7. In the other example, two sets of experiments are conducted on hybrid flow traffics generated by Skype, MSN, and GTalk with results shown in FIGS. 8 and 9. For each set of the experiments, Receiver Operating Characteristics (ROC) curves are used as the performance metric. ROC curve is a curve of detection probability, P_D, vs. false probability, P_FA, where,
P_DiΔP({circumflex over ()}_i=1|_i=1),  (38)
P_FAiΔP({circumflex over ()}_i=1|_i=0),  (39)
for ∀i=0, 1, 2. By tuning the parameters θ_σ, θ_A, and θ_V according to equation (5) and the function AUDIOVIDEOClassify described above ROC curves may be generated as shown in FIGS. 6-9.

[0063] During the experiments, network traffic flows are collected from three applications, i.e., Skype, MSN, and GTalk. For each application, network traffic flows are collected in two scenarios. For the first scenario, two laboratory computers acting as peers were communicating with each other. There was direct connection between the two peers. For the second one, a firewall was used to block direct connection between the two peers so that the application was forced to use relay nodes. For example, real time classification method is performed on first 10 seconds of each collected net work flow, i.e., _max≦10 seconds. As an example the sampling interval is set as T_s=0.5 milliseconds.

[0064] In one example, the scenario where each Skype flow carries one type of traffic is considered. In other words, for each network traffic flow, none or one of ₀, ₁, and ₂ is 1. FIG. 6A-6C show the ROC curve of classifying FILE, AUDIO, and VIDEO traffic flows.

[0065] In another example, hybrid Skype traffic flows is considered. In other words, for each network traffic flow flow, none or one or several of ₀, ₁, and ₂ may be 1. FIG. 7A-7C show the ROC curves of classifying FILE, AUDIO, and VIDEO flows, respectively.

[0066] In yet another example, experiments are performed on network traffic generated by Skype, MSN, and GTalk, as these are very common applications at the current time. In other words, an AUDIO traffic flow may be a VoIP flow generated by Skype, MSN, or GTalk. So are VIDEO or FILE traffic flows. Note that, GTalk does not support video conference. Similarly, two sets of experiments are conducted, one on single-typed network traffic flows and the other on hybrid flows.

[0067] In one example, the scenario where each flow carries one type of traffic is considered. By tuning the thresholds, θ_σ, θ_A, and θ_V, ROC curves of classifying FILE, AUDIO, and VIDEO traffic flows may be generated. FIG. 8A-8C show the ROC curves of classifying FILE, AUDIO, and VIDEO flows in single typed traffic flows from multiple network operators in accordance with one embodiment of the present invention.

[0068] In still another example, experiments on hybrid flows are performed. FIG. 9A-9C show the ROC curves of classifying FILE, AUDIO, and VIDEO flows in hybrid flows from multiple network operators in accordance with one embodiment of the present invention.

[0069] To better demonstrate FIGS. 6-9, some typical values of P_D and P_FA pairs are shown in Table II below.

[0070]

TABLE II
Typical PD and PFA values.
	Skype	Skype + MSN + GTalk
PFA (PD)	Single	Hybrid	Single	Hybrid
FILE	.039(.989)	.044(.977)	.041(.981)	.039(.967)
AUDIO	0(1)	0(1)	0(.995)	.002(.986)
VIDEO	0(.993)	0(.965)	0(.952)	0(.948)

[0071] From Table II, one notes that classification of AUDIO traffic is more accurate than that of VIDEO. Specifically, 100% accurate classification of Skype AUDIO flows may be achieved. This is due to the fact that AUDIO traffic has higher regularity than VIDEO does, as shown in FIGS. 10A and 10B, which illustrates the estimated PSD of two independent AUDIO traffic flows and two independent VIDEO traffic flows. One can immediately tell the dominant periodic component at 33 Hz in the AUDIO traffic flows. This frequency corresponds to the 30-millisecond IPD of the employed voice coding.

[0072] On the other hand, VIDEO PSDs have peaks at 0. It means that non-periodic component dominates in VIDEO traffic flows. One can see that PSDs of the two VIDEO traffic flows are close to each other. That is the reason that the real time classification methods described above may achieve high classification accuracy by using PSD features.

[0073] From Table II, one can tell that classification of single-typed flows is more accurate than that of hybrid flows. Mixing multiple types of traffic together may be equivalent to increasing noise. Hence, it is not surprising that classification accuracy is reduced.

[0074] One further notes that classification of Skype flows is more accurate than that of flows generated from general applications, i.e., Skype, MSN, and GTalk. Empirically, it may be found that Skype traffic flows are similar to traffic GTalk flows, but quite different from MSN traffic flows. For example, both Skype and GTalk AUDIO traffic flows have 33-millisecond IAT, whereas MSN AUDIO traffic flow has 25-millisecond IAT. When these flows are mixed together, classification accuracy is reduced. But the accuracy reduction is acceptable. Specifically, for hybrid AUDIO traffic flows at P_FA≈0. P_D is reduced from 1 to 0.986; while for hybrid VIDEO traffic, P_D is reduced from 0.965 to 0.948. This demonstrates that the above described classification method is robust. The robustness results from the fact that the subspace identification module composes multiple subspaces in the original high-dimensional feature space. As a result, PSD feature vectors of Skype and GTalk are likely to be within different subspaces than those of MSN. Therefore, it is still with high probability to make correct classification.

[0075] Advantages of the present invention may include one or more of the following. Real-time network traffic flows carrying multiple types of sessions may be classified for the purpose of QoS provision, which is of high interest to network carriers and operators. Network data packet features may be combined with packet arrival time in a novel way to establish a stochastic process model for representing a network flow. Spectral analysis techniques may be applied to this stochastic process model to extract regularities residing in real-time network traffic. Subspace decomposition followed by subspace bases identification provides a more general method for low dimensional subspace identification comparing to prior art method such as principal components analysis. Normalized distance to subspace used as the similarity metric for net work traffic flow classification provides a more general and robust method conventional nearest-neighbor and K-nearest-neighbor methods.

[0076] While the invention has been described with respect to a limited number of embodiments, those skilled in the art will appreciate that other embodiments can be devised which do not depart from the scope of the invention as disclosed herein. Accordingly, the scope of the invention should be limited only by the attached claims.