article: Our acoustical environment abounds with repetitive sounds, some of which are related to pitch perception. It is still unknown how the auditory system, in processing these sounds, relates a physical stimulus and its percept. Since, in mammals, all auditory stimuli are conveyed into the nervous system through the auditory nerve (AN) fibers, a model should explain the perception of pitch as a function of this particular input. However, pitch perception is invariant to certain features of the physical stimulus. For example, a missing fundamental stimulus with resolved or unresolved harmonics, or a low and high-level amplitude stimulus with the same spectral content–these all give rise to the same percept of pitch. In contrast, the AN representations for these different stimuli are not invariant to these effects. In fact, due to saturation and non-linearity of both cochlear and inner hair cells responses, these differences are enhanced by the AN fibers. Thus there is a difficulty in explaining how pitch percept arises from the activity of the AN fibers. We introduce a novel approach for extracting pitch cues from the AN population activity for a given arbitrary stimulus. The method is based on a technique known as sparse coding (SC). It is the representation of pitch cues by a few spatiotemporal atoms (templates) from among a large set of possible ones (a dictionary). The amount of activity of each atom is represented by a non-zero coefficient, analogous to an active neuron. Such a technique has been successfully applied to other modalities, particularly vision. The model is composed of a cochlear model, an SC processing unit, and a harmonic sieve. We show that the model copes with different pitch phenomena: extracting resolved and non-resolved harmonics, missing fundamental pitches, stimuli with both high and low amplitudes, iterated rippled noises, and recorded musical instruments.
  The proposed model consists of three main parts: (i) a cochlear model that translates auditory stimulus, sin (t), into the AN population activity, SAN (t, fCF) (Fig 1A); (ii) a sparse coding (SC) unit that represents the AN population response as a sparse (few non-zeros) set of coefficients, h. (Fig 1B); (iii) a readout unit that translates the active coefficients in h into a probability density function pdf (fp) of pitches. To compare the model with known psychoacoustic phenomena, we set the estimated stimulus’ pitch to be the maximum value of that pdf function, sp ∈ R1 (Fig 1C). This scalar represents the most plausible pitch for the particular given input signal, sin (t). Each part is briefly described in the following sections.
  To demonstrate the advantage of using sparse coding algorithms, we compared the performance of the algorithm (Eq 2) for sparse (λ = 0. 01) and non-sparse solutions (λ = 0, Least squares). The resulting vectors, hLS and hSC, for the two cases are shown in Fig 5 for the two aforementioned 30 dB SPL stimuli of sin, 2 (t) and sin, 3 (t). In the case of sin, 2 (t), there is little difference between the two solutions (Fig 5A). For the case of sin, 3 (t), which includes relatively higher (non-resolved) frequency components, the difference is more substantial (Fig 5B). Still, although the LS solutions (hLS) yield more nonzero coefficients than those of the sparse ones (hSC), the solutions can visibly be related to the harmonic structures of the two input stimuli (compare the resulting h vectors with the FTs in Fig 2B). The benefits of the sparse representation are evident when we introduce stimuli with high volume levels. Fig 6 compares the processing of a missing fundamental harmonic series with six harmonic components (k ∈ [3,8]), sin, L (t) =gL⋅ (∑k∈[3,8]sin (2π⋅f0k⋅t) ). (6) The fundamental frequency is set to f0 = 225 Hz, and the amplitude gL corresponds to either 30 dB (Fig 6A) or 90 dB SPL. In these simulations, to account for the nonlinearities of the cochlea in the case of the higher sound level, we used the more updated model of Zilany et al. [27–29]. Fig 6A shows that for the 30 dB SPL signal, the AN population response is limited to a small region around the CFs of the stimulus’ spectral components, as expected. On the other hand, for the 90 dB SPL case (Fig 6B), the AN population response is dramatically different as it spreads out along the whole cochlea. Moreover, since the response is heavily saturated, the spatiotemporal patterns of peaks and troughs for each of the AN fibers are biased relative to the case of the moderate sound levels. Despite this significant difference in the AN responses between different sound levels, psychoacoustic measurements indicate that the sensation of pitch is robust to that effect [6,44]. The output of the SC model, h, is shown in Fig 6C and 6D for the two stimuli levels, 30 and 90 dB-SPL, respectively. Each panel shows both the sparse (SC, λ = 0. 01, red) and non-sparse (LS, λ = 0. 0, blue) solutions. As can be expected, the number of nonzero components in the vector h is much smaller for the SC solution when compared with that obtained by the LS algorithm. This sparseness applies to both sound levels, but the difference between the two solutions is much more noticeable for higher sound level (compare the blue line in Fig 6C with that of Fig 6D). Specifically, the LS solutions for the two stimuli levels are fundamentally different, and there is no apparent preservation of the spectral components of the input stimulus. Thus, the LS solution is variant to sound level, as opposed to what we would have expected from a representation of a pitch in the auditory system [45]. In comparison, the SC solutions do manage to preserve their overall structure. While the two SC solutions are not identical, both have only a few non-zero terms that directly relate to the frequency components of the input stimulus. Consequently, it seems that the sparse requirement in Eq 2, at the expense of the accuracy of the ordinary LS solution, contributes to the invariant representation of the stimulus in the vector h, regardless of its sound level. In order to understand the effect of these different representations on pitch perception, we compare the resulting probability density functions (the pdf (fp) ) for the two stimulus levels and both LS and SC solutions (Fig 6E and 6F). For the low sound level (Fig 6E and 6F), the two solutions are alike, and there is no apparent benefit to using one over the other: both curves (blue and red) have peaks at the same frequencies, and the maximum probability point equals that of the stimulus' fundamental frequency, i. e. , f^0=225Hz. The result is substantially different for the 90 dB SPL amplitude. In this case, the pdf that corresponds to the solution of the LS algorithm has lost all resemblance with the stimulus' frequency components—it is just a flat, noisy curve. In comparison, the pdf that corresponds to the SC solution still has a clear indication of the original stimulus properties. In particular, the pdf curve (red line in Fig 6F) peaks at the harmonics of the fundamental f0, with a maximum peak at f^0=225Hz. Unlike, for example, the Fourier transform, the SC transform enables the use of various dictionaries that can be set according to some desired specifications. A standard option is to train a dictionary according to desired optimal constraints [35], but this is certainly not a prerequisite. For instance, when mapping patches of images into their parsimonious representations, one can choose to set the atoms of the dictionary as the basis of the discrete cosine transform, which is a straightforward and efficient choice (see for example Ch. 12 in [25]). In the current study, we checked two families of dictionaries: the first (Dsine) contains atoms created by sine stimuli, and the second (Dstack) contains atoms created from harmonics tones stimuli (harmonic stack). Specifically, each of the atoms in Dsine was produced by stimuli of one tone with random (uniformly distributed) amplitudes, and the atoms of Dstack were created by complex tone stimuli (harmonics 1st-6th) of the same moderate amplitude level (see Methods) The two dictionaries were checked with both stimuli of 45dB SPL (Fig 7A–7C) and high amplitude levels of 90dB SPL (Fig 7B–7D). All simulations had the same spectral structure as given by Eq 6, i. e. , all signals were complex tones with the 3–8 harmonics (the missing fundamental case). The results show the maximum peaks (blue dots) of the resulting pdfs that are taken from an interval of an octave around the fundamental frequency, as is the practice in psychoacoustic measurements (see for example [41,42]). Comparing the four subplots of Fig 7, the dissimilarities between the results of the two dictionaries demonstrate quantitatively different results but qualitatively similar performances, at least in the lower frequency regions. For example, the two dictionaries yield relatively good estimations of the pitches up to about 4k Hz of the first stimulus’ harmonic (the third harmonic). For frequencies higher than 4k Hz, there is more variance around the expected fundamental frequency (thick black dashed middle line). Specifically, the 90dB SPL estimations for both Dsine and Dstack seem to be horizontally spread (Fig 7B–7D). From our experience, these deviations can be reduced by using dictionaries of higher resolution (i. e. , with more cochlear channels and more atoms). However, due to the naïve structure of this model, there is a computational limit to the dictionary size that we can use. This restriction will hopefully be alleviated by a future model (see Discussion). Fig 8 shows in more detail four selected examples of f0 = 606. 4 Hz. The two h vectors in Fig 8A are taken from the moderate and high amplitudes levels of the Dsine dictionary case. The performance of the model for the 90dB level is a bit degraded compared to the 45dB level, as expected (see also Fig 6). Specifically, in the 90 dB SPL, the model estimates lower coefficients in h due to the nonlinearity of the cochlea. The resulting pdf, i. e. , the probability for a particular pitch given these SC coefficients, is shown in Fig 8B. Fig 8C and 8D show the h vectors and pdfs for the second dictionary, the Dstack. Comparing Fig 8A with Fig 8C for the respected amplitude levels shows that the two SC vectors have different coefficients. This difference is due to the particular structure of the atoms in each of the two dictionaries above. Still, in spite of this structural differences, the h vectors have nonzero coefficients over indices that represent harmonics of f0. Consequently, the resulting pdfs have qualitatively similar results and the maximum peaks appear at the same location. Thus, the result is that the SC model predicts the same pitches for these two cases. Interestingly, the harmonic sieve is designed, in principle, to be optimal in the case of Dsine (and under certain assumptions, see [19,40]). Even so, it can still be used for the Dstack dictionary and yield good results. Note that it might be that for the Dstack dictionary there is a better (in the aforementioned optimal sense) representation for the harmonic sieve, but we did not pursue this path any further. In summary, this section emphasizes the notion that choosing a dictionary can improve or reduce the performance of the model in different aspects. Thus, one emerging interesting question from the above discussion is which dictionary can be acquired in a biologically compelling manner to match psychoacoustic and physiological measurements best? In this paper, however, we do not address this issue, but we have chosen to focus on demonstrating that parsimonious representation of an auditory information can explain relatively high cognitive tasks, i. e. , the percept of a pitch. In what follows we assume a dictionary that is built of a single tone. Hopefully, the choosing of such a relatively simple dictionary, instead of a more intricate one, would prove to be clear and emphasize the qualitative abilities of such an approach. The harmonics of a periodic signal are spatially distributed along the BM. Because of BM properties, low harmonics create separate peaks that are translated into distinct excitation patterns in the activities of the ANs. Since higher harmonics, on the other hand, do not yield such distinct peaks, these harmonics do not have distinct excitation patterns. Consequently, low harmonics are referred to as resolved and higher harmonics, approximately at the 5th–10th harmonics [14], as unresolved. For unresolved stimuli, the temporal aspects of AN response convey more information about the pitch than the spatial aspects of that response. Thus, using stimuli with resolved and unresolved harmonics is a controlled way to inspect the temporal processing of a pitch in the auditory system. Broadly speaking, temporal models, such as the summary autocorrelation function (SACF) [4], disregard the resolvability of the stimulus. This is because such models account for the interactions between the harmonic components of the periodic signal whether they are resolved or not [44]. This indifference, however, stands in contrast to psychoacoustic observations [11,12,14,46]. The current model combines temporal and spatial aspects of the AN response within its atoms. Thus, it is interesting to examine the model’s response to this class of stimuli. We compare five stimuli of complex tones (Fig 9), sin (t) =g45dB⋅ (∑k=rr+3sin (2πf0k⋅t) ). (7) In Eq 7, each stimulus has the same spectral structure, that is, four consecutive harmonics, but the spectral locations of harmonics vary. Specifically, the spectral location of the first harmonic in each signal is set by r ∈ {1,6, 10,17,22}. Additionally, the fundamental frequency is configured to f0 = 443 Hz, and the gains (g45 dB) of all stimuli are equivalent to 45 dB SPL. Fig 9A shows the sparse coefficient vectors h for all five stimuli. Additionally, we compare the SCs of the two dictionaries, Dsine (lines) and Dstack (dashed lines). Dsine consists of tone-atoms and Dstack consists of complex tones that contain six harmonics with decreasing amplitudes (1 to 1/6). All the atoms are created with 60 dB SPL stimuli. As can be seen, there is a slight difference between the results, but overall both dictionaries yield the same SCs, pdfs, and pitch estimations. For that reason, we focused the rest of the analysis on the Dsine results. In Eq 7, Low values of r represent stimuli with resolved harmonics while larger values (r > 5) represents stimuli with unresolved harmonics. Observe (Fig 9A) that for most of the stimuli, the four prominent spectral components of h are successfully reconstructed. Still, there is an apparent degradation (from resolved to unresolved stimuli) in the amplitude of the coefficient terms. For resolved stimuli, e. g. , the 1st–4th harmonics stimulus (r = 1), the h vector holds prominent terms equivalent to the frequency components of the input signal. As r increases, the estimated terms in h decrease; for example, in the 22nd–25th harmonics stimulus (r = 22), the vector h only contains two out of the four equivalent frequency components of the input signal (Fig 9A, inset). This degradation is due to the reduced ability of the AN population to phase-lock with high frequency stimuli. Thus, the match between the atoms and the stimuli is less accurate, which results in smaller SC coefficients. Fig 9B shows the corresponding pdfs of each of the five stimuli. The right panel shows the entire pdf while the left panel focuses on the vicinity of the fundamental f0 = 443 Hz. The numbers in the figure indicate the arrangement in a descending order of the peak heights, starting from the highest peak (peak number 1). For all stimuli, the maximum of each density function was obtained at the fundamental frequency. It is instructive to note that for each of the stimuli, other local maxima, indicating other pitch possibilities, are allocated at harmonic ratios of f0. For example, in the case of the lowest harmonic complex (r = 1), maximum peaks in the pdf are also available at the f0/4, f0/2, and 2f0. These other options represent the octave equivalence of the perceived pitch. Usually, humans perceive these options to be the same pitch, or to have the same pitch chroma. Hence, in psychoacoustic measurements, this mixing between octaves is generally not considered as an error [41–43] (specifically, see SI in [41]). In Fig 10 we compare the pdfs of the two complex tones of Fig 9. The stimulus that contains the 1–4 harmonics is shown in blue while the stimulus that contains the 22–25 harmonics is in green. This comparison is done over one octave to avoid the other pitch equivalence solutions (which are approximately at the same height of the 1st peak). By inspection (Fig 10A), the ratio between the 1st and the 2nd peaks is higher for the resolved stimulus (blue line) compared to that of the unresolved stimulus (green line). We chose to denote this difference in the ratios as the salience of the stimuli. Fig 10B shows a comparison of additional stimuli with different fundamental frequencies; each one is a complex tone that contains four consecutive tones. The colors of the circles match the stimuli colors of Fig 9. Note that the salience of all stimuli decline with the increase of the location of the first harmonic, i. e. , as the harmonics transcend from resolved to unresolved. Finally, stimuli with high harmonics also have additional peaks around the locus of the harmonic components (see the fLOCUS peaks in Fig 9B). These phenomena are consistent with known physiological data [9]. To conclude, despite the model’s seemingly spatially-based nature, it can derive the pitch of unresolved harmonics to some extent. Unlike purely temporal models, however, it penalizes these stimuli in relation to resolved ones. Note that this penalty is a consequence of the cochlear properties and not of the SC module. Specifically, this penalty was not introduced artificially into the model—it is an implicit property of the atoms and stems directly from the properties of the cochlear model. Perceived pitches are usually considered within the context of periodic signals. For example, the perceived pitch of a complex tone is its fundamental frequency, f0, whether it exists in the complex or not. Consider the following harmonic series sin (t) =g⋅ (∑k∈[k0, k1]sin (2π (f0k+Δf) ⋅t) ), (8) For Δf = 0, human subjects usually perceive the pitch of sin (t) as f0, the fundamental frequency of the harmonic signal [47]. This is true even for cases when f0 is not present in the signal, i. e. , k0 > 1. For increased Δf > 0, the expression in Eq 8 is no longer harmonic. Nonetheless, previous psychoacoustic experiments have revealed that human subjects do manage to perceive pitches with these shifted stimuli, and the detected pitches are approximately shifted on a linear scale relative to the fundamental frequency [9,40,47–49]. These kind of stimuli are important because they demonstrate that pitch detection does not follow the stimulus' envelope, which does not change in this case (this is not true for stimuli with unresolved harmonics); nor does it follow the spacing between frequency components of a stimulus [1]. Thus, the phenomenon of pitch shift was used as a counter example for models that exploit the temporal envelope of a stimulus or other of its temporal features, such as zero crossing, peaks, etc. [18]. Fig 11A–11D shows the model' s solutions of the sparse coefficient vectors, h, for different shifted signals (Eq 8). In these simulations, the input signal, sin (t), has four frequency components that are set to k ∈ [4,7], and the (missing) fundamental frequency is set to f0 = 200 Hz. In this example, setting Δf = 0 creates a stimulus sin (t) that is a complex harmonic series of f0 = 200 Hz. For a shift of Δf = 40 Hz, the stimulus is no longer a complex tone of f0 = 200 Hz, and a shift of Δf = 100 Hz changes the stimulus to be a complex tone of f0 = 100 Hz (with the [9,11,13,15] harmonics). Finally, for a frequency shift of Δf = 200 Hz, the input signal is once again a complex tone of the fundamental f0 = 200 Hz, but this time with the [5,6, 7,8] harmonics. Following these observations, we would expect the model to exhibit this ambiguity by the frequency shift Δf, and to alternate its predictions between the frequencies of f0,0. 5f0,2f0, etc. We performed 500 simulations of sin (t) (Eq 8) with f0 = 200 Hz. Each signal contains the first six harmonics (i. e. , k ∈ [1,6]), and each is simulated with a different Δf. Fig 11E shows the estimated frequency f^0 as a function of the lower harmonic component fL in the input signal. For each of the stimuli, the corresponding pdfs are calculated, and the four highest peaks are indicated. We chose to include the four prominent peaks of each pdf to show the ambiguities of these signals. The estimated pitches are clustered in lines as a function of the lower frequency fL, consistent with the known “first effect of pitch shift” [47]. One can also see that as the shift Δf increases, the slopes of the estimated lines slightly decrease, in accordance with the known psychoacoustic phenomenon [47]. To conclude, the current model qualitatively reproduces the known psychoacoustic phenomenon of pitch shift, even though these aperiodic signals are not part of the model' s dictionary. This implies that the current model: generalizes to complex new stimuli; it does not depend on the stimulus temporal envelope for cues; and it does not use the spacing between the stimulus harmonics to estimate pitch. Transposed tones (TTs) were first introduced to explore the relative sensitivity of the auditory system for binaural timing stimuli [50]. Oxenham et al. [13] used these signals to check the relation between the tonotopic organization in the cochlea and the perception of pitch. The motivation of using these signals is to introduce low-frequency temporal structures into the basal part of the cochlea that usually processes high frequencies. TTs are produced by the modulation of half-wave rectified sine waves with carrier waves. Due to the limited synchronization and low-pass properties of the basilar membrane, the outer hair cells, and the ANs, the fine details of the carrier waves would be negligible. Thus, the results are half-wave rectified sine waves of low frequency (f0) in high CFs regions (> 4 kHz). In this manner, there is a separation between spatial locations and fine temporal structures along the cochlea. Measurements carried out by Oxenham et al. [13] showed that subjects could not estimate the fundamental frequency of the TT stimuli. This inability means that the spatial arrangement of CFs along the cochlea is essential for the perception of pitch. Oxenham et al. have also shown that the summary autocorrelation function (SACF), a well-known temporal model, is indifferent to TTs and thus to the tonotopic organization in the cochlea. When applying TTs to the analysis of the SACF, the model does manage to extract the correct fundamental frequencies, in contrast to the aforementioned psychoacoustic evidence. In the case of the proposed model, each region along the cochlea is characterized by its local spatiotemporal activities. These localized patterns are embedded in the atoms for each CF. Hence, we predicted that the SC model would not be able to separate between spatial and temporal processing. Fig 12 shows the processing of TT stimuli by the SC model. We first considered three TTs with the fundamentals f0 = 230 Hz, 250 Hz, and 270 Hz. The expression for each of these TTs is given by sTT (t) =g30dB⋅∑k∈[1,3]s¯k (t) ⋅sin (2πfc, k⋅t). (9) In this equation, s¯k (t) is a low-pass filtered version of a rectified sine wave (see Methods), and each of these tones is modulated by the three carriers: fc, 1 = 4 kHz, fc, 2 = 6. 35 kHz, and fc, 3 = 10. 08 kHz. Fig 12A shows the results of the sparse coefficient vectors (h) for the above three TT stimuli. It is apparent that the three sparse representations occupy the same indices, i. e. , the same frequencies fd in h. Therefore, the SC model cannot distinguish between these stimuli based on their sparse representations. And indeed, the pdfs of these sparse vectors are inconclusive (Fig 12B). Next, we simulated a batch of 100 epochs (Fig 12C). Each epoch contained the TT of Eq 9, and each had an incremented fundamental frequency taken from the interval f0 ∈ [100 Hz, 500 Hz] (red bars in Fig 12C). We repeated the simulation also for pure tones (blue bars in Fig 12C) and compared the two by normalizing the measurements with the respective f0. The results are consistent with the findings of Oxenham et al. , that is, the SC model could not estimate the f0 for the TT stimuli successfully. Delaying a signal of broadband white noise and adding it back to the original one creates a signal known as rippled noise. When this process of delaying and adding is repeated, a signal known as iterated rippled noise (IRN) is created. These signals contain temporal regularities in the time domain and spectral peaks at the reciprocal of the delayed time in the spectrum domain. Due to the nature of these signals, human listeners report perceiving two sensations: a tonal part that amounts to the pitch of the reciprocal of the delay (d ms) and an additional noisy sensation [51,52]. Repeating the iteration process results in a more prominent sensation of the tonal pitch [53]. Adding the delayed noisy signal back to the original one with a gain of one (delay-add) yields a signal with spectral peaks that are located at the reciprocal of the delay time d. But adding the delayed signal with a gain of minus one (delay-subtract) yields a signal with peaks in the power spectrum that are shifted by 1/2d, as if the delay-subtract signal is an odd-harmonic complex of half the frequency of the delay-add version. Delay-add stimulus raises a sensation of pitches of 1/d Hz, whereas delay-subtract is usually perceived to be more ambiguous and yields pitches that are slightly higher or lower than 1/d Hz [53]. Simulations of delay-add and delay-subtract stimuli with the delays of d = 2,4, and 5 ms, are shown in Fig 13. These signals are created as follows: for a white noise, x (t), the iterated signal sn (t) is created by si (t) =si−1 (t) +g⋅si−1 (t−d) s0 (t) =x (t), (10) for i = 1, …, nitr, and nitr is the number of iterations (e. g. , 1,2, or 10). All simulations were done using Carney' s model (Zilany et al. [27–29]) with stimuli levels of 70dB SPL; the dictionary contained 1000 sine-atom groups and each group has 10 time-shifts (phases) in it (g = 10, Eq 3). In the first row (Fig 13A–13C), which contains the delay-add cases, a clear peak appears at the reciprocal of the delay, 1/d, as expected. As the number of repetitions increases, so does the prediction quality of the model, i. e. , more estimations are concentrated around 1/d Hz. In Fig 13D–13F we show the delay-subtract cases for the same delays. In these simulations (Fig 13E and 13F), the inferred pitches are located around the reciprocal of the delays, as expected. But we would also expect the measurements to peak at approximately 1/d±10%, which does not happen. Music in the Western culture is based on a musical scale that relates periodic (or quasiperiodic) sounds to their fundamental frequencies. Thus, musical instruments that are based on this musical scale produce harmonic sounds based on these fundamental frequencies. As such, different instruments have different spectral coloring (i. e. , timbre), but human listeners can perceive and compare the fundamental frequencies between the instruments [54]. This ability is due to the pitch perception property of clustering periodic (or quasiperiodic) sounds into classes, i. e. , musical notes. In this section, we checked the SC model with recorded musical notes [55]. For this, we used a dictionary with 1000 atoms, each of which had ten different phases (g = 10, Eq 3). Each recorded stimulus was divided into Tsteps = 100 time steps that were analyzed separately. Fig 14A shows the FT of a recorded violin note of A5 (880 Hz) played with a bow (arco). At each time step Tsteps, Eq 3 is solved separately to obtain h˜∈RgM. Next, the coefficients of each group are summed together to get the SC vector h ∈ RM. The matrix Hg∈RM×Tsteps is the aggregation of all h ∈ RM over Tsteps = 100 steps (Fig 14B). Each of the columns of Hg (the SC vectors) are then processed by the harmonic sieve unit to produce the probability of that time step. The collection of all these pdf is given by the matrix Pg∈RM×Tsteps (Fig 14C). Finally, to have one single probability for each stimulus, the matrix Pg is averaged over the time domain and normalized appropriately (Fig 14D). As in previous cases, the pitch of the signal is defined as the maximum point in this pdf. We repeated this procedure with recorded notes of a flute, a violin, and a piano (Fig 15). All results are shown on a chromatic scale. Each dot in Fig 15 is the estimated pitch of a recorded instrument; the colored text indicates the played note.
  The SC model presented here combines both spatial and temporal aspects of the AN population response. On the one hand, the SC model is based on atoms that have limited spatial support, namely, the nonzero section along the BM that is given by the equivalent rectangular bandwidth of the cochlea [60]. On the other hand, each atom also includes temporal information about the activities of the AN fibers at that spatial location. The current model is structured in a similar manner to pattern-matching models [5,16,18]. The Fourier-like spectrum analyzer that extracts the resolved harmonic components of a given stimulus [18] is modified to include both the cochlea and the SC modules. Both models use templates to associate an estimated pitch with presented stimuli [58,61]. Despite these similarities, there is a fundamental difference between pattern-matching models and the current one. The atoms (templates) used in this study contain both spatial and temporal activities (Fig 3). Consequently, although the patterns of the AN population activities may be spectrally unresolved, there is still enough spatiotemporal information for identifying the different harmonics (Fig 9). Additionally, the SC does not rely on synthetic bases such as the sine and cosine of the Fourier transform (FT) but actual AN activities. Of course, these bases are also unresolved for high tones which means that the model exhibits less salience with unresolved stimuli, again, in line with psychoacoustic experiments [13]. Note that since the SC model implicitly inherits this property for high tones from the known attributes of the cochlear model, there is no need to add this feature into the model explicitly. Historically, models of pitches were linked to harmonic analysis theories. Thus, it is instructive to note the mathematical connection between the SC model and the FT. Indeed, for the particular case of deprecating the dictionary matrix D into a square matrix, and setting λ = 0, the optimal solution for any periodic signal x of Eq 2 is given by the FT. In this case, the optimal solution is for h to be the coefficients of the FT and D its matrix [62]. In the SC case, applying λ ≠ 0 enables the use of biologically oriented, non-orthogonal, and redundant dictionaries [63]. In this sense, the current proposed model can be seen as an extension of classic pattern-matching models. How does the SC model compare with temporal ones? Traditional temporal models, the most prominent of which are based on the summary autocorrelation function (SACF), also exploit temporal features from the AN population responses [2,3, 10,64]. However, there are several open issues with these type of models. First, SACF models need to have long tapped delay lines for the correlation module, i. e. , about 40ms and maybe more if noise is accounted for [15,65]. Currently, however, there is no physiological evidence to support such structures [9,31]. In contrast, the SC model exploits local spatiotemporal features without the need for long tapped delay lines. For that reason and to keep the model biologically plausible, we chose to use only short time segments of 5ms (see [66]). Different temporal interval durations were also tested. For shorter time intervals, e. g. , 2ms, the SC model acts as a place-rate model; i. e. , it managed to estimate only low resolved signals by their spatial activities along the BM. Longer time intervals also improved the model predictions for unresolved stimuli, up to a maximal estimation of about 10ms. Another issue associated with temporal models is that they treat resolved and unresolved frequencies in a similar manner in contrast to known psychoacoustic measurements [11,12,14,46]. This may imply that tonotopic organization is not necessary for auditory processing, but, again, physiological evidence suggests otherwise [13]. Tonotopic organization is preserved in the auditory system up to the auditory cortex for all mammals [67–70], and current evidence suggests that pitch processing is also sensitive to it [9,14]. Recently, Laudanski et al. [71] proposed a structural theory of pitch that considers both the spatial and the temporal aspects of the AN population response. Within this framework, the perception of pitch is derived from correlated activity in pairs of points in the spatiotemporal representation of AN activities. These two points are not necessarily located along the temporal activity of the same AN fiber (pure temporal processing), nor between different AN fibers at a particular time (pure spatial processing). Both the structural theory and the proposed SC model are strongly related. In the SC model we incorporate the so-called cross-channel delays of the structural theory in the spatiotemporal patterns of the atoms. Specifically, cross-channel delays of a stimulus are compared to other cross-channel delays of the model that are embedded in the atoms. We think, however, that approaching the problem of pitch estimation from the SC aspect offers considerable benefits. First, the SC approach provides a mathematical framework that generalizes to other modalities whereas the structural theory approach offers a more specific pointwise approach [23]. Second, cross-channel delays of the structural theory can be simply acquired under the SC method by using predesigned atoms (as shown in this paper) or unsupervised training of atoms (see [36,72–74] to name just a few). These techniques were already tested, including within the auditory system [57], with great success. Thus, the SC framework can explain different possible options for such cross-channel correlations. Other theories and models that exploit the spatial, the temporal, or both, include: the spatial cross-correlation theory of Loeb et al. [75]; de Cheveigne’s solution for the problem of tapped delay lines in temporal theories [65]; Carney’s model of phase-opponency [76], Shamma et al. ’s lateral-inhibition and cross-correlation-matrix model [77–79], and the MASD of Cedolin et al. [80]. Note, however, that, these models consider only a small subset of the whole two-dimensional spatiotemporal structure created by the AN fibers. For example, Loeb et al. proposed comparing two locations along the BM that vibrate with the same phase, that is, a spatial comparison without the time domain. De Cheveigne proposed to compensate phase shifts between adjacent cochlear filters, i. e. , extracting temporal lags and discarding the spatial information; and Cedolin et al. proposed a model that is based on spatial derivation between cochlear filters with a temporal summation, namely, accounting for the differences between two adjacent cochlear filters and averaging over the time domain. Additionally, these models account well for resolved stimuli but not for unresolved ones [71]. The use of the harmonic sieve can be considered from different perspectives. First, from the probabilistic point of view it can be seen as an implementation of a likelihood function: the probability of a particular pitch given a set of (parsimonious) coefficients h. This approach originates from the pattern-matching theory, and since the proposed model can be seen as an extension to the pattern-matching models, the same theoretical and experimental motivations also apply here. For example, following Goldstein et al. [16], we chose the templates of the harmonic sieve to be Gaussian functions [16,40,58,81]. It might be that for different dictionaries, e. g. , dictionaries that contain harmonic stacks, there are better options, but we were not concern with optimizing this feature in this paper. Second, from the physiological perspective, the harmonic sieve can be thought of as a simple feedforward neural network. In such instance, a set of Gaussian templates of one tone (one row in the matrix G) can be seen as a neuron with a modulated selectivity curve, i. e. , a neuron that responds to a particular tone and its successive octaves. For examples of such implementations, see [61,82]. Finally, the harmonic sieve can be considered as a simple (i. e. , linear) readout function that extracts the perceived pitch from the activities of the spatiotemporal receptive fields and introduces it in a manner that enables an easy comparison with psychoacoustic data. Third, from the biological perspective, it had been shown [83] that harmonic templates of this sort can emerge naturally from basic processing in the auditory periphery. Specifically, Shamma et al. demonstrated that the fundamental features include: frequency analysis, fast changing delays at the CFs, phase-locking, and half-wave rectification. All of these properties part of the cochlear models that we used. In this paper, we checked two types of dictionary families: the first is created by sine stimuli while the second is created by stimuli of harmonic tones (i. e. , harmonic stacks). These dictionary types were constructed and tested for various amplitude levels. Choosing a dictionary can influence the congruity of the SC model’s results with psychoacoustic measurements: the sensitivity of the model to resolved and unresolved stimuli, the response to low and high stimuli levels, etc. All these features emerge from the cochlear properties that are encapsulated within the atoms. Thus, an important question is which dictionary can be acquired in a biologically compelling manner and will best match psychoacoustic measurements? We intend to investigate this interesting question further in a future paper. However, in this paper we focused on the main premise—that parsimonious representation of auditory information can explain relatively high cognitive tasks, such as that of the percept of a pitch. Accordingly, we chose to work, for the most part, with a simple dictionary of pure tones. Indeed, pure sinusoidal stimuli are rare in our natural acoustic surroundings [84] and there is no guarantee that the auditory system has access to such components at all. However, note that there is a subtle difference between plausible stimuli for a (hypothetical) training process and the outcome of this process, the atoms themselves. This is an important distinction because it implies that it is reasonable to assume that a learning process over natural-like stimuli, for example, vowels and consonants, can yield local spatiotemporal atoms and not necessarily stack-like atoms. Moreover, there is circumstantial physiological evidence to support similar spatiotemporal structures along different areas in the auditory system. For example, Norman-Haignere et al. [85] researched specific regions in the anterior auditory cortex that strongly react to resolved harmonic tones, and, to a lesser extent, to unresolved ones. Additionally, Carney et al. found, in the anteroventral cochlear nucleus [66], cells that have distinct spatiotemporal tuning patterns in response to pure tones. Acquiring a dictionary directly from the AN population response is not the necessarily the only implementation nor is it the optimal one. A different approach that has been successfully applied and has many variants [24,34–37] is to perform an unsupervised training from a randomly selected stimulus set that represents sampled statistics of the input domain. These unsupervised learning algorithms create dictionaries that are oriented to specific kind of inputs, for example, natural sounds [73]. Subsequently, they usually generate dictionaries that lead to sparser and more accurate results. It is thus important to realize that the proposed SC technique is not limited to a particular set of atoms (see for example Ch. 12 in [25]). Throughout this paper, we tried to keep the discussion at the representational level[59]. Accordingly, we did not introduce a plausible neural network mechanism to concentrate on what we saw as the central theme of the current paper—the generalized principle of using sparse representation also for high perception tasks—such as the estimation of pitches. We felt sufficiently confident to follow such a path because the current literature already includes several plausible neural network implementations of sparse coding [86–90] Another important point is that we did not introduce a state-of-the-art solution but a qualitative one. For this reason, the current model operates on stimuli with fixed time intervals. One relatively simple and standard technique to alleviate this restriction is by running the same model in consecutive times (for example, [73]). In the current paper, we normalized all AN population activities. By normalization we mean that we divided each AN simulated result by its maximum response. The normalization that we propose stems from the assumption that relevant information about the pitch is related to the overall spatiotemporal structures of the AN population responses and not their absolute instantaneous rate level. The problem of keeping the estimation of the pitch invariant to the stimulus’ level is due to other deformations and nonlinearities in the AN responses: saturation in the activity of the AN fibers; change in the locations of peaks of ANs population activities; or the relative phases between the different AN fibers [31,45]. It is important to stress that this normalization neither changes nor corrects these effects. It is worthwhile noting that normalization has been observed across the central neural system in general and in the auditory system in particular [91,92]. This stems from the fact the different modalities need to process a large dynamic range of stimuli, whether it is brightness in the visual sensors or a change of few orders of magnitude in the level of the sound stimulus.
  We showed that sparse coding principles that were successfully applied to other modalities can explain pitch perception. This general approach of a parsimonious representation of the sensory information is the main premise of this paper and this finding resonates with ideas of a canonical computation by the nervous system [23,74,93–98]. Specifically, sparse representation of information can explain neural activity in the visual cortex [24,34,95,99]; the olfactory system of insects [100,101]; and findings in the mammalian auditory cortex [102–105] (however, see [106]). Hopefully, this paper is one small step in searching for such a generalized theory.
  We used two types of cochlear models: (i) Slaney' s MATLAB auditory toolbox [26] and (ii) the more moderate model of Zilany et al. [27–29]. We used Slaney’s model for most of the simulations in this paper because it is relatively accurate and computationally fast for moderate sound levels. This model, however, does not account for the high amplitude level nonlinearities of the cochlea; for this phenomenon we used Zilany’s et al. model (Figs 5 and 6). In Zilany’s model, the parameters were set as: (1) the outer and inner hair cells were taken to be in healthy condition; (2) the fractional Gaussian noise that is related to the spike rate generated by each AN was approximated in order to save computation time; (3) the model was originally built to match the AN population response of the cat, but we used the built-in option to tune it to the human cochlea [60]; and (4), only high spontaneous rate AN fibers were used. For both cochlear models, the input stimulus and the simulation of the AN fibers were performed with a sampling rate of FS = 100k Hz; both had the same number of AN fibers (N = 200). The simulations of the AN population responses were performed as follows: a 15ms stimulus, sin (t) ∈RTa, is constructed and the AN population response is calculated using the chosen model. Finally, the number of samples is Ta = 15ms∙FS (rounded to an integer if necessary). The output of the cochlear model is a matrix of time-samples over the number of AN fibers, i. e. , Ta × N. This matrix is then truncated to contain only the last 5ms samples, i. e. , SAN (t, fCF) ∈ RT×N and T = 5ms∙FS. Lastly, SAN (t, fCF) is normalized by its maximum value (0 < |SAN (t, fCF) |<1). Since the model of Zilany et al. has upper and lower bounds on possible simulated CFs, between 125Hz to 20k Hz, all simulations were performed within this frequency interval. Finally, following the usual SPL format, the stimulus levels are introduced into the cochlea in pascal units normalized by the threshold of hearing, i. e. , gPa=20μ⋅10gdB/20. Each of the atoms dj ∈ RT⋅N×1, j ∈ [1, M], in D are taken as the AN population response to a pure tone. Each dj is the vector form of the respective AN population response matrix, i. e. , the vectorized form. The vectorization is performed using MATLAB’s convention of stacking the matrix’s columns one after the other. The dictionary matrix D = [d1, …, dM] ∈ RT⋅N×gM contains M = 1000 atoms for each group g. The groups are collection of atoms with the same CFs but with different phases (Eq 3). Note that D is a rectangular (N∙T ≫ M), highly redundant, matrix. The number of atoms was set arbitrarily by trial-and-error. From our experience, fewer atoms (e. g. , M = 250), also yielded reasonable results. In this paper we chose to implement the SC (Eq 2) by means of the least absolute shrinkage and selection (LASSO) algorithm [32]. We did so because it has a simple implementation, has a relatively acceptable running time, and usually yields good results. The implementation of the LASSO involves an iterative solution derived by gradient descent. Specifically, the vector h is the solution of the following iterated equation: hk+1=soft (hk+1αDT (vAN−Dhk), λ2α). (11) In this equation, vAN is the vectorized form of the AN population response SAN (t, fCF) (after vectorization and normalization), and the operator soft (x, T) = sign (x) max{0, |x| − T} is defined for each entry in the vector x. The algorithm runs until a convergent criterion is met or until a pre-set number of iterations is exceeded. For the algorithm to converge, the parameter α should maintain a certain condition (α ≥ max eig (DTD) ). For the simulations in this paper we used the LASSO implementation within MATLAB Inc. [107]. We also created a tweaked version of this algorithm, but there was no substantial difference between the two implementations. The pitch estimation unit is a variant of the known harmonic sieve [58] implementation. It denotes the likelihood of a particular pitch given the sparse coefficient vector h ∈ RM. We implemented it as a multiplication between h˜∈RP and the matrix G ∈ RP×P (Eq 4) in which P = 15M, interpolating from the M values of h. Each row in G corresponds to a candidate pitch fp ∈ [125Hz, 20kHz] and is composed of a set of Gaussian weights at successive harmonics for this particular pitch, i. e. , Gi, j=∑kexp (− (fj−nkfi) 2σi22). (12) The index i is taken from 1 to ⌊ (20k−125) /σp⌋. Thus, if the sparse vector h has a coefficient that relates to the harmonics of the fi pitch, multiplying by the matrix G would emphasize (give high score) to that entry in the pdf. Otherwise, if h does not contain harmonics that relate to the fi pitch, the result would be a low score in the pdf. The standard deviations of each of the Gaussian curves in G is a function of fi, σi=0. 2⋅ERB (fi) =0. 2⋅24. 7 (4. 37fi+1). (13) In this equation, ERB stands for the equivalent rectangular bandwidth of Glasberg and Moore [60]. We tried different variations of σi, including piecewise curves, all with relatively similar qualitative results. The main constraint was to avoid overlap between the Gaussian distributions and to keep adjacent Gaussian curves wide enough to account for noise in the pitch cues in h˜. Finally, using this scheme we only had to calculate the matrix G once, making the algorithm relatively efficient and fast. The dictionary Dsine ∈ RT⋅N×M is constructed with stimuli of one tone. It has N = 300 CF channels and M = 1500 atoms, each formed by a tone stimuli of uniformly selected amplitudes over the interval of 30dB to 70dB SPL (gsine∼unif (30dB, 70dB) ), ssine (t) =gsine⋅sin (2π⋅f0k⋅t) (14) The dictionary Dstack has 200 CF channels and 1000 atoms. Each of these atoms contains a stimulus of six consecutive harmonics (1st-6th); these harmonics have linearly decreasing amplitudes (from 1 to 1/6), sstack (t) =∑k=16gstackk⋅sin (2π⋅f0k⋅t). All of the Dstack atoms are formed by gstack = 40dB SPL level stimuli. All simulations were made using Carney' s model (Zilany' s et al. [27–29]) to account for these high amplitude levels. In order to avoid aliasing in the spatiotemporal domain, the maximum frequency in Dstack is the maximum frequency (20k Hz) divided by the number of harmonics (6). The transposed tone stimuli in this paper are modulated rectified tones with three carrier waves: fc, 1 = 4k Hz, fc, 2 = 6. 35k Hz, and fc, 3 = 10. 08k Hz (Eq 9). The modulated tone is given by s¯k (t) =[sin (2πf0⋅t) ]+*hLP, k (t). (15) In this equation, the rectification operator is given as [x]+ = max (0, x). The operation * is a convolution, and hLP, k (t) is a four-order Butterworth low-pass filter (see Oxenham et al. paper [13] for more details). The cutoff frequencies of the low-pass filter is taken as 0. 2fc, k for each of the three modulated frequencies (k = {1,2, 3}). Following Yost et al. [42], all the stimuli have amplitude levels of 70dB SPL. To account for this amplitude level, all simulations were performed using Carney' s model (Zilany et al. [27–29]). To account for the random phases of the stimuli we used g = 10 groups (Eq 3). The same simulations with a dictionary that contains no groups, g = 1 and ϕg = 0, have been slightly noisier but qualitatively the same. Each of the six cases shown in Fig 13 is a histogram of 500 simulations. All the repetition cases are normalized respectively. The maximum peaks of the pdfs are selected within an interval of one octave around the pitch frequency, 1/d. All simulations are performed using the same dictionary; this dictionary contains atoms that are 5ms long and have one tone at a level of 45dB SPL. Finally, all stimuli are filtered by a low-pass with a frequency band of 4k Hz [42]. We used a dictionary with 1000 sine-atoms of length 5ms. Each such sine-atom composed a group of 10 different phases (g = 10, Eq 3). For these simulations, we used Slaney' s MATLAB toolbox [108] as the cochlear model (much faster); the cochlea had 300 CF channels. We used Eq 2 for each of the Tsteps = 100 timesteps separately. Each of the SC coefficient vectors h˜∈RgM was averaged and normalized appropriately into h ∈ RM (each group separately). The collection of all these vectors, for all Tsteps, formed the matrix Hg∈RM×Tsteps. Each of these SC vectors (the columns of Hg) are then processed by the harmonic sieve to produce the probability of pitch Pg∈RM×Tsteps at each time step (Fig 14C). Finally, to have one single probability for each stimulus, the matrix Pg is averaged over the time domain and normalized appropriately (Fig 14D). As in previous cases, the pitch of the signal is defined as the maximum point in this pdf. All measurements were downloaded from the University of Iowa, electronic music studio, from the musical instrument samples page [55].
summary: By means of a sound' s pitch, we can easily discern between low and high musical notes, regardless of whether they originate from a guitar, piano or a vocalist. The relation between different sounds that yield the same percept is what makes pitch an interesting subject of research. Today, despite extensive research, the mechanism behind this physical to perceptual transformation is still unclear. The large dynamic range of the cochlea combined with its nonlinear nature makes the modeling and understanding of this process a challenging task. Given a large amount of physiological and psychological data, a general explanation consistent with many of these phenomena would be a major step in elucidating the nature of pitch perception. In this paper, we recast the problem in the general framework of sparse coding of sensory stimuli. This framework, initially developed for the visual modality, posits that the goal of the neural representation is to represent the flow of sensory information in a concise and parsimonious way. We show that applying this principle to the problem of pitch perception can explain many perceptual phenomena.
keywords: medicine and health sciences
  music cognition
  ears
  psychoacoustics
  social sciences
  neuroscience
  hearing
  simulation and modeling
  perception
  inner ear
  cognitive psychology
  research and analysis methods
  sensory physiology
  head
  auditory system
  psychology
  anatomy
  cochlea
  physiology
  psychophysics
  biology and life sciences
  sensory systems
  sensory perception
  cognitive science
  music perception
  pitch perception
title: A New Approach to Model Pitch Perception Using Sparse Coding