Personalized HRTF generation method based on spectrum feature extraction and clustering analysis

Abstract

The invention discloses a personalized HRTF generation method based on spectrum feature extraction and clustering analysis, and belongs to the technical field of three-dimensional audio, the method comprises two parts of feature extraction and feature clustering analysis on an HRTF spectrum, and the method comprises the following steps: firstly, loading and pre-processing input sofa data, hilbert transform and a cross-correlation function are used to acquire a binaural time difference ITD, an amplitude difference is calculated to acquire a binaural sound intensity difference ILD, sub-band division of an equivalent matrix bandwidth is used to determine a center frequency, direction transfer function data DTF is acquired through a designed Gammatone filter, positive gradient information is calculated, and finally, a relation matrix and a clustering center are acquired through features, so that the clustering accuracy is improved. And carrying out dimension reduction mapping by using MDS, and obtaining a clustering result and a clustering center by using an algorithm to realize three-dimensional sound field rendering. The method is not limited by the audiometry environment of the user, the personalized HRTF of the user is autonomously, conveniently and quickly matched, accurate positioning is achieved, confusion is reduced, and the method can be popularized and used in scenes such as virtual reality and games.

Description

Technical Field

[0001] The present invention relates to the field of three-dimensional audio technology, and in particular to a personalized HRTF generation method based on spectrum feature extraction and cluster analysis. Background Art

[0002] Head Related Transfer Function (HRTF) simulates the complex acoustic effects of diffraction, reflection, and scattering of sound waves caused by the physiological structures of the head, auricle, and trunk when sound waves propagate from the sound source to the human ears. Therefore, HRTF can be used to simulate the state of sound source signals in different spatial orientations when they are transmitted to the human ears, thereby simulating a realistic binaural three-dimensional sound effect.

[0003] Commonly used methods for personalized HRTF include acoustic measurement method, simulation calculation method, selection method, regression method, and perceptual feedback method.

[0004] The personalized HRTF method based on acoustic measurement needs to obtain IR (Impulse Response) in an anechoic room or semi-anechoic room environment, which has high requirements on the environment and professionalism, and the measurement process is complicated and time-consuming.

[0005] Simulation calculation is to obtain HRTF by calculating the sound pressure law of the sound source sound wave at the entrance of the human ear canal through numerical calculation to simulate the sound wave appearing as the boundary of the free field on the human body surface. This method requires the use of 3D scanning and other equipment, which is difficult to implement and difficult to popularize commercially.

[0006] The selection method is based on matching the user's human morphological features with the morphological features in the existing HRTF sample library, and finding the most similar features as the user's personalized HRTF. However, it cannot guarantee the best subjective listening and speaking effect when rendering the audition;

[0007] The regression method generates personalized HRTF by establishing a relationship mapping between human morphological features and HRTF. Since there is no unified standard for the representation capabilities of different human morphological features, some signals in HRTF cannot be extracted from human body measurement features or ear images. The data in the HRTF dataset is limited, and the measurement standards of each dataset are not unified, so they cannot be combined and used.

[0008] The perceptual feedback method is to test a large number of HRTF samples through human ears, and the user subjectively selects the HRTF sample with the best test effect as the user's personalized HRTF. However, according to the physiological characteristics of human hearing, long-term testing of a large number of samples can easily lead to auditory fatigue, which can cause confusion or unclear distinction of the later test audio, thus affecting the effectiveness of HRTF.

[0009] To summarize the currently commonly used personalized HRTF methods, the above methods have problems such as harsh measurement environment, poor subjective audiometric effect during feature matching, inaccurate audio positioning due to long-term audio measurement fatigue, and the generation of confusing sounds. Summary of the invention

[0010] The present disclosure provides a personalized HRTF generation method based on spectrum feature extraction and cluster analysis, so that the user can achieve the best effect of subjective audio-visual three-dimensional sound effect through independent short-term audio measurement without being restricted by the audio measurement environment.

[0011] A method including HRTF spectrum feature extraction and cluster analysis is proposed, wherein the HRTF spectrum feature extraction step A includes the following steps:

[0012] Step A1: load the input sofa data to obtain sofa_obj, and normalize the data in sofa_obj;

[0013] Step A2: passing the data sofa_obj obtained in step A1 into a low-pass filter, and using Hilbert transform and cross-correlation function to obtain the spectral feature interaural time difference itd;

[0014] Step A3: for the data sofa_obj obtained in step A1, obtain the frequency spectrum feature binaural sound intensity difference ild by calculating the amplitude difference;

[0015] Step A4: Divide the data sofa_obj obtained in step A1 into sub-bands to determine the center frequency fc;

[0016] Step A5: Design and use a Gammatone filter for sofa_obj obtained in step A1 and fc obtained in step A4 to obtain a directional transfer function dtf;

[0017] Step A6: Calculate the positive gradient information extraction of the directional transfer function dtf obtained in step A5 to obtain the positive gradient data pge and the center frequency gfc of the positive gradient data;

[0018] A personalized HRTF generation method based on spectrum feature extraction and cluster analysis is modeled, including spectrum feature extraction and cluster analysis of HRTF, wherein the cluster analysis step B of HRTF includes the following steps:

[0019] Step B1: passing the feature parameters obtained in the HRTF feature extraction module (step A) to the calM function to obtain the relationship matrix M and the cluster center W;

[0020] Step B2: Use the MDS (Multidimensional Scaling) algorithm to perform dimensionality reduction mapping on the relationship matrix M obtained in step B1 to obtain M';

[0021] Step B3: Use the K-Means algorithm to cluster M' obtained in step B2, obtain the clustering index idxs, and obtain the clustering result and cluster center through idxs;

[0022] Among them, the normalization in step A1 is to normalize both the radius and HRTF data, wherein the radius normalization is to set the radius value to 1, and the HRTF normalization is to use the maximum value of the impulse response when the subject's azimuth and elevation are 0° for normalization; the low-pass filter uses a Butterworth filter, which obtains the maximum auditory cross-correlation by using the filter, and then uses the Hilbert transform and cross-correlation function to obtain the interaural time difference.

[0023] This technology provides a personalized HRTF generation method based on spectrum feature extraction and cluster analysis. Compared with acoustic measurement and simulation calculation, this method requires less listening time and is not limited by the environment and equipment. Compared with the selection method and regression method, this technology allows users to select the best subjective listening effect through listening. In addition, in order to solve the problem of auditory fatigue caused by excessive listening,

[0024] Feature extraction is performed through the existing HRTF, and then the HRTF is divided into several categories using the clustering method. The user only needs to listen to the HRTF at the center of the audition category, and then listen to the HRTF in the selected category (feature extraction and clustering can also be performed again in the category for subjective listening). Compared with a single audition, the number of HRTFs that the user needs to listen to is very small, which reduces the time required for the user to select the best HRTF during the subjective listening process. Due to the larger perceived difference in the cluster center, the user's listening effect is significant.

[0025] In this solution, users can reduce the HRTF rendered audio required for perception, reduce auditory fatigue, and large differences can easily select similar HRTF data. BRIEF DESCRIPTION OF THE DRAWINGS

[0026] Figure 1 is a flowchart of a personalized HRTF generation method based on spectrum feature extraction and cluster analysis according to the present disclosure;

[0027] Figure 2 It is a schematic diagram of a user interface UI for performing a perception test according to an embodiment of the present disclosure. DETAILED DESCRIPTION

[0028] Now we will refer to the following in detail: Figure 1The figure shows a flow chart of a personalized HRTF generation method based on spectrum feature extraction and cluster analysis of the present invention, which mainly includes HRTF spectrum feature extraction and feature cluster analysis. The main steps include: loading and preprocessing the input sofa data, using Hilbert transform and cross-correlation function to obtain binaural time difference, calculating amplitude difference to obtain binaural sound intensity difference, sub-band division of equivalent matrix bandwidth to determine center frequency, design and use of Gammatone filter to obtain directional transfer function, calculate positive gradient information, obtain relationship matrix and cluster center through features, use MDS for dimensionality reduction mapping, and use K-Means algorithm to obtain clustering results and cluster centers. The specific implementation details of each step are as follows:

[0029] Step A1: First load the data:

[0030] sofa_obj = SOFAload (sofa)

[0031] SOFAload() is the function for loading HRTF data, and sofa is the HRTF file in the data set.

[0032] Then normalize:

[0033] features.coord.pos=sofa_obj.SourcePosition

[0034] features.coords.pos(:,3)=1

[0035] idx=find(ismember(feature.coords.pos,[azi ele r],'rows'))

[0036] sofa_frontal=sofa_obj.Data.IR(idx,:,:)

[0037] features.IR=sofa_obj.Data,IR. / (max(abs(sofa_frontal(:)))+eps)

[0038] Among them, SourcePosition is an attribute containing azimuth, elevation, and radius. Normalizing the radius is to directly set the radius value to 1; ismember() is to determine whether the row (azi ele r) in the array is included. If the row is (azi ele r), the row value is 1 when returned, otherwise it is 0; find() is used to find the index of non-zero elements in the array. The code uses ismember() and find() functions to find the index of the specified azimuth and elevation. This step is to find the index of the front position, that is, the position where the azimuth and elevation are both 0 (the radius is normalized, so the radius is 1); sofa_obj.Data.IR is the impulse response (IR) of HRTF, abs() is an absolute value operation, max() obtains the maximum value in IR, eps is a minimum value in MATLAB to prevent the divisor from being 0. Here, all IRs are normalized by obtaining the maximum value in the IR in front and performing dot division.

[0039] Step A2: Obtain the interaural time difference itd from the data sofa_obj obtained in step A1.

[0040] [lp_a,lp_b]=butter(butterpoly,cut_off_freq_norm)

[0041] f_ir=filter(lp_a,lp_b,squeeze(sofa_obj.Data.IR(ii,jj,:)))

[0042] e_sir1=abs(hilbert(squeeze(f_ir(ii,1,:))))

[0043] e_sir2=abs(hilbert(squeeze(f_ir(ii,2,:))))

[0044] [IACC(ii),idx_lag]=max(abs(xcorr(e_sir1,e_sir2)))

[0045] itd(ii)=idx_lag-sofa_obj.API.N

[0046] features.itd=itd / sofa_obj.Data.SamplingRate

[0047] IR is 1550×2×256 data, 1550 is the number of spatial orientations, 2 is left and right ears, and 256 is the IR data of a single ear in each orientation, so ii in the above is processing one orientation at a time, and jj is processing one ear at a time; sofa_obj.Data.SamplingRate is the sampling rate; butter() is the constructed Butterworth low-pass filter, where butterploy is the order, cut_off_freq_norm is the standardized cutoff frequency, and the return values ​​are the polynomial coefficients of the numerator and denominator respectively; filter() is to use the constructed filter to filter the signal; squeeze() simplifies the dimension of the data; hilbert() is the Hilbert transform function to obtain the envelope information; xcorr() is to calculate the time difference between the two signals, where the signals are left and right ear signals, and the time difference between the left and right ears is calculated; finally, dividing by the sampling rate is to convert the number of samples into seconds and save them in features (Note: all the features obtained below will be saved in features, which may not be written in some descriptions).

[0048] Step A3: Convert the amplitude values ​​of the left and right ears into decibel values ​​dB, and obtain the binaural sound intensity difference ild by calculating the difference between the left and right ears.

[0049] ild=mag2db(left)-mag2db(right)

[0050] Where left and right are the amplitudes of the HRTF, and the amplitudes are converted to decibels using the mag2db() function.

[0051] Step A4: Calculate the equivalent matrix bandwidth (ERB) scale,

[0052] audlimits=9.2645*sign(freq).*log(1+abs(freq)*0.00437)

[0053] audrange=audlimits(2)-audlimits(1)

[0054] n=floor(audrange / bw)

[0055] audpoints=audlimits(1)+(0:n)*bw+(audrange-n*bw) / 2

[0056] fc=(1 / 0.00437)*sign(audpoints).*(eps(abs(audpoints) / 9.2645)-1)

[0057] Where freq is a 2×1 array containing the lower and upper frequency limits ([700,18000] is used in the implementation example), audlimits is the ERB scale corresponding to freq, converted from Hertz to the equivalent matrix bandwidth (ERB) scale; audrange is the number of ERBs in the freq range; n is the number of ERB subbands divided in the calculation, bw is the number of ERBs counted as one subband (the implementation example uses bw=1), and the floor() function rounds down; audrange-n*bw is the remaining ERB after the calculation of the subband division, audpoints is the center frequency point after the frequency range freq is divided according to the equivalent matrix bandwidth (ERB) scale; fc is the conversion from the ERB scale back to frequency.

[0058] Step A5:

[0059] dtf' = 2*real(Gammatone(IR))

[0060] dtf”=sqrt(max(squeeze(permute(reshape(dtf',dim1),dim2)),0))

[0061] dtf=20*log10(norm(squeeze(dtf”)) / sqrt(size(dtf”)(4)))

[0062] Where Gammatone() is the designed 4th order Gammatone filter,

[0063] Gammation() is the designed 4th-order gammatone filter, which filters the impulse response IR of HRTF data; real() is the operation to solve the real part of the complex number; dim1 and dim2 are the dimensions of the transformed shape; reshape() reshapes the tensor, norm() performs normalization operation; sqrt() takes the square root; 20*log(·) converts the result into decibels (dB).

[0064] Step A6:

[0065] pge(ii,jj)=max(dtf(ii,jj+1,:)-dtf(ii,jj,:),0)

[0066] gfc=sqrt(fc(1:size(dtf,2)-1).*fc(2:end))

[0067] Among them, ii, jj are the sizes of the first dimension and the second dimension of dtf respectively, and each value is processed one by one; fc is the center frequency in step A4; size() obtains the size of a certain dimension.

[0068] All features in step A will be saved in features and transferred to step B for feature clustering.

[0069] Step B1: The specific implementation of the calM() function is shown below

[0070] [IDXs,C,sumd,D]=kmeans(feature,N,'Replicates',R)

[0071] M=setM(string(A(IDXs==ii)),A,M,alpha)

[0072] [m,W]=min(D)

[0073] Among them, the parameters accepted by the calM() function are feature, M, N, A, R, alpha; features is the feature obtained in step A, the features are: itd, ild, dtf, IR, pge, gfc and other required information, feature is one of these features, one of these features is passed into the function for clustering, after clustering according to a certain feature, clustering is performed according to other features; M is the initialized relationship matrix, which is used to store the relationship strength value between two objects; N is the number of categories to be classified; A is the name of all objects; R is the number of clustering rounds; alpha is the weight coefficient; k The means() method is used for clustering, and the returned results are: IDXs is the cluster label of each data point, C is the location of the cluster center, sumd is the sum of the distances from each data point to the center of its category, and D is the distance from each data point to the centers of all categories; the setM() method calculates the relationship matrix after clustering, saves the result in M ​​and returns it, where ii is the value of 0-N in each loop, and the loop is used to calculate the relationship for each category; the min() method returns the minimum value and index of each category, because the cluster center distance is 0, so the location of the cluster center is obtained; the calM() function returns the relationship matrix M and the cluster center index W.

[0074] Step B2:

[0075] MM=ones(size(M,1),size(M,2)).*ma

[0076] d=MM-M

[0077] n=size(d,1)

[0078] t(i,j)=-0.5*(d(i,j)^2-1 / n*d(i,:)'-1 / n*d(:,j)'*d(:,j)+i / n^2*sum(sum(d.^2)))

[0079] [V,D]=eig(t)

[0080] XX=V(:,1:2)*D(1:2,1:2).^(1 / 2)

[0081] M'=[real(XX(:,2)),imag(XX(:,1))]

[0082] Among them, size() is used to obtain the size of a certain dimension; ma is a constant 4; d is to map the relationship matrix M to d, where the stronger the relationship, the closer the corresponding object distance; t is the centralized distance matrix, which can represent the relative distance relationship between objects. The centralization process eliminates the deviation of the data and obtains a matrix suitable for dimensionality reduction processing; i, j are traversed in the double-layer loop for each value; eig() calculates the eigenvalue and eigenvector of the matrix, V is the corresponding eigenvector, and D is the eigenvalue; M' is the coordinate value of each object mapped to the low dimension.

[0083] Step B3:

[0084] [idxs,C,sumd,D]=kmeans(M',N,'Replicates',R)

[0085] [m,W]=min(D)

[0086] Among them, idxs stores the category corresponding to the index, and W is the index of the category center.

[0087] This technology provides a personalized HRTF generation method based on spectrum feature extraction and cluster analysis. Compared with acoustic measurement and simulation calculation, this method takes less time, and the equipment only needs headphones and players, and does not involve any professional operation. Compared with the selection method and regression method, the user's subjective listening and speaking is best selected directly through perception. The selection method and regression method cannot guarantee the best results in subjective listening and speaking, and users are more confident in their own listening results. Similarly, in order to solve the problem of auditory fatigue caused by excessive user audition, the existing HRTF is extracted by using features, and then the HRTF is divided into several categories by clustering. The user only needs to audition the HRTF at the center of the category, and then audition the HRTF in the selected category (feature extraction and clustering can also be performed again in the category for subjective listening). Compared with a single audition, the number of HRTFs that the user needs to listen to is very small, and the user can quickly select the HRTF with the best perception. Moreover, the difference in perception of the cluster center is large, and the user is more likely to hear the difference. In this solution, the user can reduce a lot of HRTF rendering audio that needs to be perceived, reduce auditory fatigue, and the large difference can easily select similar HRTF data.

[0088] Although embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions and variations may be made to the embodiments without departing from the principles and spirit of the present invention, and that the scope of the present invention is defined by the appended claims and their equivalents.

Claims

1. A personalized HRTF generation method based on spectrum feature extraction and cluster analysis, characterized in that: The personalized HRTF generation method based on spectrum feature extraction and cluster analysis includes: The spectral feature extraction step A of HRTF includes the following steps A1 to A6: Step A1: Load the input sofa data to obtain sofa_obj, and normalize the data in sofa_obj; Step A2: input the data sofa_obj obtained in step A1 into a low-pass filter, and use Hilbert transform and cross-correlation function to obtain the interaural time difference itd; Step A3: for the data sofa_obj obtained in step A1, use the amplitude difference calculation method to obtain the binaural sound intensity difference ild; Step A4: Divide the data sofa_obj obtained in step A1 into sub-bands to determine the center frequency fc; Step A5: design a Gammatone filter for the sofa_obj obtained in step A1 and the center frequency fc obtained in step A4 to obtain a directional transfer function dtf; Step A6: extracting the positive gradient of the directional transfer function dtf to obtain the positive gradient data pge and the center frequency gfc of the positive gradient data; The cluster analysis step B of HRTF includes the following steps B1 to B3: Step B1: The characteristic parameters obtained in step A are input into the calM function to obtain the relationship matrix M and the cluster center W; Step B2: Based on the multidimensional scaling algorithm MDS, the relationship matrix M obtained in step B1 is subjected to dimensionality reduction mapping to obtain a mapping matrix M'; Step B3: Based on the K-Means algorithm, cluster the mapping matrix M' obtained in step B2, obtain the clustered index idxs, and obtain the clustering result and cluster center through the index idxs; Among them, the normalization in step A1 is to normalize both the radius and HRTF data, wherein the radius normalization is to set the radius value to 1, and the HRTF normalization is to use the maximum value of the impulse response when the subject's azimuth and elevation are 0° for normalization; the low-pass filter in step A2 is a Butterworth filter, which uses the filter to obtain the maximum auditory cross-correlation, and then uses the Hilbert transform and cross-correlation function to obtain the interaural time difference.

2. The personalized HRTF generation method based on spectrum feature extraction and cluster analysis according to claim 1, characterized in that: The step A1 comprises the following steps: Data loading is based on the following formula: sofa_obj = SOFAload (sofa) Among them, sofa_obj is HRTF data, SOFAload() is the function of loading HRTF data, and sofa is the HRTF file in the data set. Normalization is performed based on the following formula: features.coord.pos=sofa_obj.SourcePosition features.coords.pos(:,3)=1 idx=find(ismember(feature.coords.pos,[azi ele r],'rows')) sofa_frontal=sofa_obj.Data.IR(idx,:,:) features.IR=sofa_obj.Data,IR. / (max(abs(sofa_frontal(:)))+eps) Among them, SourcePosition is an attribute that includes azimuth, elevation, and radius. Normalizing the radius is to directly set the radius value to 1; ismember() is used to determine whether the row azi ele r in the array is included. If the row is azi ele r, the row value is 1 when returned, otherwise it is 0; find() is used to find the index of non-zero elements in the array. The code uses ismember() and find() functions to find the index of the specified azimuth and elevation. Here, the index of the position directly in front is found, that is, the radius is normalized at the position where the azimuth and elevation are both 0, so the radius is 1; sofa_obj.Data.IR is the impulse response (IR) of HRTF, abs() is an absolute value operation, math() obtains the maximum value in IR, eps is a minimum value in MATLAB to prevent the divisor from being 0. Here, all IRs are normalized by obtaining the maximum value in the IR directly in front and performing dot division.

3. The personalized HRTF generation method based on spectrum feature extraction and cluster analysis according to claim 1, characterized in that: The step A2 comprises the following steps: Based on the following formula, itd is obtained for the data sofa_obj obtained in step A1: [lp_a,lp_b]=butter(butterpoly,cut_off_freq_norm) f_ir=filter(lp_a,lp_b,squeeze(sofa_obj.Data.IR(ii,jj,:))) e_sir1=abs(hilbert(squeeze(f_ir(ii,1,:)))) e_sir2=abs(hilbert(squeeze(f_ir(ii,2,:)))) [IACC(ii),idx_lag]=max(abs(xcorr(e_sir1,e_sir2))) itd(ii)=idx_lag-sofa_obj.API.N features.itd=itd / sofa_obj.Data.SamplingRate Among them, IR is 1550×2×256 data, 1550 is the number of spatial orientations, 2 is left and right ears, and 256 is the IR data of a single ear in each orientation, so ii in the above is processing each orientation one by one, and jj is processing each ear one by one; sofa_obj.Data.SamplingRate is the sampling rate; butter() is the constructed Butterworth filter, where butterploy is the order, cut_off_freq_norm is the standardized cutoff frequency, and the return values ​​are the polynomial coefficients of the numerator and denominator respectively; filter() is to use the constructed filter to filter the signal; squeeze() simplifies the dimension of the data; hilbert() is the Hilbert transform function to obtain the envelope information; xcorr() is to calculate the time difference between the two signals, where the signals are left and right ear signals, and the time difference between the left and right ears is calculated; finally, dividing by the sampling rate is to convert the number of samples into seconds and save it in features.

4. The personalized HRTF generation method based on spectrum feature extraction and cluster analysis according to claim 1, characterized in that: The step A3 comprises the following steps: The amplitude values ​​of the left and right ears are converted into decibel values ​​dB, and the binaural sound intensity difference ild is obtained by calculating the difference between the left and right ears. The calculation method is as follows: ild=mag2db(left)-mag2db(right) Where left and right are the amplitudes of the HRTF, and the amplitudes are converted to decibels using the mag2db() function.

5. The personalized HRTF generation method based on spectrum feature extraction and cluster analysis according to claim 1, characterized in that: The step A4 comprises the following steps: audlimits=9.2645*sign(freq).*log(1+abs(freq)*0.00437) audrange=audlimits(2)-audlimits(1) n=floor(audrange / bw) audpoints=audlimits(1)+(0:n)*bw+(audrange-n*bw) / 2 fc=(1 / 0.00437)*sign(audpoints).*(eps(abs(audpoints) / 9.2645)-1) Among them, freq is a 2×1 array containing the lower and upper frequency limits, audlimits is the ERB scale corresponding to freq, converted from Hertz to the equivalent matrix bandwidth ERB scale; audrange is the number of ERBs in the freq range; n is the number of ERB subbands divided in the calculation, bw expresses the number of ERBs required to calculate one subband, and the floor() function rounds down; audrange-n*bw is the remaining ERB after calculating the subband division, audpoints is the center frequency point after the frequency range freq is divided according to the equivalent matrix bandwidth ERB scale; fc is converted back from the ERB scale to frequency.

6. The personalized HRTF generation method based on spectrum feature extraction and cluster analysis according to claim 1, characterized in that: Step A5 includes the following steps: dtf1 = 2 * real (Gammatone (IR)) dtf2=sqrt(max(squeeze(permute(reshape(dtf',dim1),dim2)),0)) dtf=20*log10(norm(squeeze(dtf”)) / sqrt(size(dtf”)(4))) Gammation() is a pre-designed 4th-order gammatone filter that filters the impulse response IR of the HRTF data; real() is an operation to solve the real part of the complex number; dim1 and dim2 are the dimensions of the transformed shape; reshape() reshapes the tensor, norm() performs normalization; sqrt() takes the square root; and 20*log(·) converts the result into decibels (dB).

7. The personalized HRTF generation method based on spectrum feature extraction and cluster analysis according to claim 1, characterized in that: The step A6 comprises the following steps: pge(ii,jj)=max(dtf(ii,jj+1,:)-dtf(ii,jj,:),0) gfc=sqrt(fc(1:size(dtf,2)-1).*fc(2:end)) Among them, ii, jj are the sizes of the first and second dimensions of dtf respectively, and each value is processed one by one; fc is the center frequency in step A4; size() obtains the size of a certain dimension. All features in step A will be saved to features for feature clustering in step B.

8. The personalized HRTF generation method based on spectrum feature extraction and cluster analysis according to claim 1, characterized in that: The step B1 comprises the following steps: The specific implementation of the calM() function is shown below: [IDXs,C,sumd,D]=kmeans(feature,N,'Replicates',R) M=setM(string(A(IDXs==ii)),A,M,alpha) [m,W]=min(D) Among them, the parameters accepted by the calM() function are feature, M, N, A, R, alpha; feature is the feature passed in by the calling function, the feature is obtained by step A, the features are itd, ild, dtf, pge, gfc, the features in this information are passed to the kmeans() function for clustering, clustering according to a certain feature and then clustering according to other features; M is the initialized relationship matrix, which is used to store the relationship strength value between two objects; N is the number of categories to be classified; A is the name of all objects; R is the number of clustering rounds; alpha is the weight coefficient; kmeans( ) method is used for clustering, and the returned results are: IDXs is the cluster label of each data point, C is the location of the cluster center, sumd is the sum of the distances from each data point to the center of its category, and D is the distance from each data point to the center of all categories; the setM() method calculates the relationship matrix after clustering, saves the result in M ​​and returns it, where ii is the value of 0-N in each loop, and the loop is used to calculate the relationship for each category; the min() method returns the minimum value and index of each category, because the cluster center distance is 0, so the location of the cluster center is obtained; the calM() function returns the relationship matrix M and the cluster center index W.

9. The personalized HRTF generation method based on spectrum feature extraction and cluster analysis according to claim 1, characterized in that: The step B2 comprises the following steps: MM=ones(size(M,1),size(M,2)).*ma d=MM-M n=size(d,1) t(i,j)=-0.5*(d(i,j)^2-1 / n*d(i,:)'-1 / n*d(:,j)'*d(:,j)+i / n^2*sum(sum(d.^2)))[V,D]=eig(t) XX=V(:,1:2)*D(1:2,1:2).^(1 / 2) M'=[real(XX(:,2)),imag(XX(:,1))] Among them, size() is used to obtain the size of a certain dimension; ma is a constant 4; d is to map the relationship matrix M to d, where the stronger the relationship, the closer the corresponding object distance; t is the centralized distance matrix, which represents the relative distance relationship between objects. The centralization process eliminates the deviation of the data and obtains a matrix suitable for dimensionality reduction processing; i, j are traversed in a double-layer loop for each value; eig() calculates the eigenvalue and eigenvector of the matrix, V is the corresponding eigenvector, and D is the eigenvalue; M' is the coordinate value of each object mapped to the low dimension.

10. The personalized HRTF generation method based on spectrum feature extraction and cluster analysis according to claim 1, characterized in that: The step B3 comprises the following steps: [idxs,C,sumd,D]=kmeans(M',N,'Replicates',R) [m,W]=min(D) Among them, idxs is the category corresponding to the index, and W is the index of the category center.

Images