Rowan-Classes/6th-Semester-Spring-2024/DSP/Labs/FinalProject/obj_evaluation/comp_fwseg_mars.m
2024-04-25 18:38:09 -04:00

494 lines
19 KiB
Matlab

function [SIG,BAK,OVL]= comp_fwseg_mars(cleanFile, enhancedFile);
% ----------------------------------------------------------------------
% MARS Frequency-variant fwSNRseg objective speech quality measure
%
% This function implements the frequency-variant fwSNRseg measure based
% on MARS analysis (see Chap. 10, Sec. 10.5.4)
%
%
% Usage: [sig,bak,ovl]=comp_fwseg_mars(cleanFile.wav, enhancedFile.wav)
%
% cleanFile.wav - clean input file in .wav format
% enhancedFile - enhanced output file in .wav format
% sig - predicted rating [1-5] of speech distortion
% bak - predicted rating [1-5] of noise distortion
% ovl - predicted rating [1-5] of overall quality
%
%
% Example call: [s,b,o] =comp_fwseg_mars('sp04.wav','enhanced.wav')
%
%
% References:
% [1] Chapter 10, Sec 10.5.4,
% [2] Chapter 11
%
% Authors: Yi Hu and Philipos C. Loizou
% (critical-band filtering routines were written by Bryan Pellom & John Hansen)
%
% Copyright (c) 2006 by Philipos C. Loizou
% $Revision: 0.0 $ $Date: 10/09/2006 $
% ----------------------------------------------------------------------
if nargin~=2
fprintf('USAGE: [sig,bak,ovl]=comp_fwseg_mars(cleanFile.wav, enhancedFile.wav)\n');
fprintf('For more help, type: help comp_fwseg_mars\n\n');
return;
end
[data1, Srate1, Nbits1]= wavread(cleanFile);
[data2, Srate2, Nbits2]= wavread(enhancedFile);
if ( Srate1~= Srate2) | ( Nbits1~= Nbits2)
error( 'The two files do not match!\n');
end
len= min( length( data1), length( data2));
data1= data1( 1: len)+eps;
data2= data2( 1: len)+eps;
wss_dist_matrix= fwseg( data1, data2,Srate1);
wss_dist=mean(wss_dist_matrix);
SIG= sig_mars( wss_dist( 1), wss_dist( 2), wss_dist( 3), wss_dist( 4), ...
wss_dist( 5), wss_dist( 6), wss_dist( 7), wss_dist( 8), ...
wss_dist( 9), wss_dist( 10), wss_dist( 11), wss_dist( 12), ...
wss_dist( 13), wss_dist( 14), wss_dist( 15), wss_dist( 16), ...
wss_dist( 17), wss_dist( 18), wss_dist( 19), wss_dist( 20), ...
wss_dist( 21), wss_dist( 22), wss_dist( 23), wss_dist( 24), ...
wss_dist( 25));
SIG=max(1,SIG); SIG=min(5, SIG); % limit values to [1, 5]
BAK= bak_mars( wss_dist( 1), wss_dist( 2), wss_dist( 3), wss_dist( 4), ...
wss_dist( 5), wss_dist( 6), wss_dist( 7), wss_dist( 8), ...
wss_dist( 9), wss_dist( 10), wss_dist( 11), wss_dist( 12), ...
wss_dist( 13), wss_dist( 14), wss_dist( 15), wss_dist( 16), ...
wss_dist( 17), wss_dist( 18), wss_dist( 19), wss_dist( 20), ...
wss_dist( 21), wss_dist( 22), wss_dist( 23), wss_dist( 24), ...
wss_dist( 25));
BAK=max(1,BAK); BAK=min(5, BAK); % limit values to [1, 5]
OVL= ovl_mars( wss_dist( 1), wss_dist( 2), wss_dist( 3), wss_dist( 4), ...
wss_dist( 5), wss_dist( 6), wss_dist( 7), wss_dist( 8), ...
wss_dist( 9), wss_dist( 10), wss_dist( 11), wss_dist( 12), ...
wss_dist( 13), wss_dist( 14), wss_dist( 15), wss_dist( 16), ...
wss_dist( 17), wss_dist( 18), wss_dist( 19), wss_dist( 20), ...
wss_dist( 21), wss_dist( 22), wss_dist( 23), wss_dist( 24), ...
wss_dist( 25));
OVL=max(1,OVL); OVL=min(5, OVL); % limit values to [1, 5]
%-------------------------------------------------
function Y= bak_mars( FWSEG_VA, V5, V6, V7, V8, V9, V10, V11, V12, ...
V13, V14, V15, V16, V17, V18, V19, V20, ...
V21, V22, V23, V24, V25, V26, V27, V28)
BF1 = max(0, V21 - 0.282);
BF2 = max(0, FWSEG_VA + 9.094);
BF3 = max(0, - 9.094 - FWSEG_VA );
BF5 = max(0, 10.089 - V11 );
BF7 = max(0, 3.624 - V26 ) * BF3;
BF8 = max(0, V24 - 5.584) * BF5;
BF9 = max(0, 5.584 - V24 ) * BF5;
BF10 = max(0, V19 - 8.030) * BF1;
BF11 = max(0, 8.030 - V19 ) * BF1;
BF12 = max(0, V27 - 4.858) * BF1;
BF13 = max(0, 4.858 - V27 ) * BF1;
BF14 = max(0, FWSEG_VA + 7.282) * BF1;
BF15 = max(0, - 7.282 - FWSEG_VA ) * BF1;
BF17 = max(0, 9.458 - V16 ) * BF10;
BF18 = max(0, V27 - 10.431) * BF11;
BF19 = max(0, 10.431 - V27 ) * BF11;
BF21 = max(0, 11.059 - V22 ) * BF1;
BF22 = max(0, V26 - 8.675) * BF1;
BF23 = max(0, 8.675 - V26 ) * BF1;
BF25 = max(0, 11.195 - V6 ) * BF10;
BF26 = max(0, V8 - 7.138) * BF1;
BF27 = max(0, 7.138 - V8 ) * BF1;
BF29 = max(0, 9.006 - V10 ) * BF26;
BF30 = max(0, V14 - 8.210) * BF15;
BF35 = max(0, 7.026 - V19 ) * BF15;
BF36 = max(0, V11 - 3.424) * BF27;
BF39 = max(0, 5.418 - V17 ) * BF23;
BF40 = max(0, V28 - 6.813);
BF41 = max(0, 6.813 - V28 );
BF42 = max(0, V26 - 5.998) * BF14;
BF43 = max(0, 5.998 - V26 ) * BF14;
BF44 = max(0, V5 + 0.206) * BF41;
BF45 = max(0, - 0.206 - V5 ) * BF41;
BF46 = max(0, V22 - 7.901) * BF45;
BF49 = max(0, 7.496 - V8 ) * BF44;
BF51 = max(0, 7.904 - V11 ) * BF45;
BF52 = max(0, V26 - 10.938) * BF27;
BF54 = max(0, V9 - 4.507) * BF26;
BF56 = max(0, V28 - 0.549) * BF15;
BF57 = max(0, 0.549 - V28 ) * BF15;
BF58 = max(0, V25 - 3.252) * BF41;
BF59 = max(0, 3.252 - V25 ) * BF41;
BF60 = max(0, V23 - 7.650) * BF58;
BF61 = max(0, 7.650 - V23 ) * BF58;
BF62 = max(0, V25 - 9.931) * BF44;
BF63 = max(0, 9.931 - V25 ) * BF44;
BF64 = max(0, V25 - 4.923) * BF21;
BF65 = max(0, 4.923 - V25 ) * BF21;
BF67 = max(0, 3.746 - V28 ) * BF10;
BF68 = max(0, V11 - 5.346) * BF41;
BF69 = max(0, 5.346 - V11 ) * BF41;
BF70 = max(0, V12 - 9.026) * BF68;
BF71 = max(0, 9.026 - V12 ) * BF68;
BF73 = max(0, - 2.668 - V28 ) * BF21;
BF74 = max(0, V24 - 7.028) * BF41;
BF75 = max(0, 7.028 - V24 ) * BF41;
BF77 = max(0, - 0.224 - V6 ) * BF74;
BF78 = max(0, V5 - 3.884);
BF79 = max(0, 3.884 - V5 );
BF80 = max(0, V15 - 5.019) * BF78;
BF83 = max(0, - 1.880 - V28 ) * BF13;
BF84 = max(0, V7 - 3.067) * BF12;
BF85 = max(0, 3.067 - V7 ) * BF12;
BF87 = max(0, 5.353 - V6 );
BF88 = max(0, V13 - 3.405) * BF9;
BF89 = max(0, 3.405 - V13 ) * BF9;
BF91 = max(0, 5.599 - V13 ) * BF45;
BF92 = max(0, V15 - 9.821) * BF8;
BF94 = max(0, V14 + 2.594) * BF79;
BF97 = max(0, 8.635 - V23 ) * BF94;
BF99 = max(0, 1.332 - V24 ) * BF45;
BF100 = max(0, V7 - 0.209) * BF1;
Y = 2.751 + 0.135 * BF1 - 0.037 * BF2 + 0.328 * BF3 - 0.098 * BF5 ...
+ 0.988 * BF7 + 0.014 * BF8 - 0.034 * BF11 - 0.011 * BF12 ...
- 0.013 * BF13 - 0.002 * BF17 + 0.014 * BF18 ...
+ 0.004 * BF19 - 0.007 * BF21 - 0.017 * BF22 ...
- .895791E-03 * BF25 + 0.011 * BF26 - 0.009 * BF27 ...
- 0.007 * BF29 + 0.052 * BF30 + 0.022 * BF35 ...
- 0.002 * BF36 - 0.005 * BF39 - 0.059 * BF40 ...
- 0.050 * BF41 + 0.001 * BF42 + .743730E-03 * BF43 ...
+ 0.011 * BF44 + 0.022 * BF45 + 0.009 * BF46 ...
+ 0.004 * BF49 - 0.005 * BF51 + 0.010 * BF52 ...
- 0.001 * BF54 - 0.005 * BF56 - 0.015 * BF57 ...
- 0.032 * BF59 + 0.009 * BF60 - 0.002 * BF61 ...
- 0.009 * BF62 - 0.001 * BF63 + .819374E-03 * BF64 ...
+ 0.002 * BF65 + 0.003 * BF67 + 0.024 * BF69 ...
- 0.011 * BF70 - 0.004 * BF71 + 0.013 * BF73 ...
- 0.026 * BF74 + 0.005 * BF75 + 0.253 * BF77 ...
- 0.065 * BF78 + 0.014 * BF80 - 0.010 * BF83 ...
+ 0.001 * BF84 + 0.018 * BF85 - 0.050 * BF87 ...
- 0.002 * BF88 - 0.020 * BF89 + 0.003 * BF91 ...
- 0.043 * BF92 + .707581E-03 * BF97 - 0.015 * BF99 ...
- 0.005 * BF100;
function Y= sig_mars( FWSEG_VA, V5, V6, V7, V8, V9, V10, V11, V12, ...
V13, V14, V15, V16, V17, V18, V19, V20, ...
V21, V22, V23, V24, V25, V26, V27, V28)
BF1 = max(0, V7 - 9.535);
BF2 = max(0, 9.535 - V7 );
BF3 = max(0, V27 - 1.578);
BF5 = max(0, V6 - 5.422);
BF6 = max(0, 5.422 - V6 );
BF8 = max(0, 11.333 - V19 );
BF10 = max(0, - 6.774 - FWSEG_VA );
BF11 = max(0, V10 - 6.255) * BF8;
BF12 = max(0, 6.255 - V10 ) * BF8;
BF13 = max(0, V24 - 3.894);
BF15 = max(0, V5 - 3.884);
BF16 = max(0, 3.884 - V5 );
BF17 = max(0, V28 - 7.918);
BF18 = max(0, 7.918 - V28 );
BF19 = max(0, V13 - 6.077) * BF18;
BF20 = max(0, 6.077 - V13 ) * BF18;
BF22 = max(0, 6.614 - V20 ) * BF10;
BF23 = max(0, FWSEG_VA + 0.936) * BF8;
BF25 = max(0, V23 - 5.039);
BF26 = max(0, 5.039 - V23 );
BF28 = max(0, 9.007 - V20 ) * BF25;
BF29 = max(0, V25 - 7.582);
BF30 = max(0, 7.582 - V25 );
BF31 = max(0, V11 + 3.336) * BF16;
BF32 = max(0, V26 - 1.877);
BF35 = max(0, - 5.749 - FWSEG_VA ) * BF6;
BF36 = max(0, V7 - 4.451) * BF29;
BF37 = max(0, 4.451 - V7 ) * BF29;
BF38 = max(0, V14 - 10.158);
BF39 = max(0, 10.158 - V14 );
BF41 = max(0, 7.172 - V17 ) * BF39;
BF43 = max(0, 7.810 - V24 ) * BF26;
BF44 = max(0, V8 + 1.636) * BF3;
BF45 = max(0, FWSEG_VA - 10.068) * BF39;
BF47 = max(0, V23 - 4.721) * BF30;
BF48 = max(0, 4.721 - V23 ) * BF30;
BF50 = max(0, - 2.397 - V24 ) * BF16;
BF51 = max(0, V14 - 1.428) * BF17;
BF53 = max(0, V16 + 1.940) * BF18;
BF54 = max(0, V10 - 9.442) * BF18;
BF56 = max(0, V10 + 2.144) * BF16;
BF58 = max(0, 1.969 - V26 ) * BF2;
BF59 = max(0, V19 - 6.089) * BF16;
BF62 = max(0, 8.952 - V21 ) * BF15;
BF63 = max(0, V24 - 7.371) * BF3;
BF65 = max(0, V22 - 8.908) * BF6;
BF66 = max(0, 8.908 - V22 ) * BF6;
BF67 = max(0, V27 - 9.485) * BF30;
BF69 = max(0, V18 - 8.608) * BF10;
BF71 = max(0, V13 - 3.374) * BF25;
BF73 = max(0, V14 - 3.616) * BF13;
BF75 = max(0, V18 - 10.321) * BF32;
BF76 = max(0, 10.321 - V18 ) * BF32;
BF78 = max(0, 3.972 - V15 ) * BF26;
BF79 = max(0, V14 - 7.105) * BF26;
BF80 = max(0, 7.105 - V14 ) * BF26;
Y = 2.638 - 0.089 * BF1 + 0.083 * BF5 - 0.162 * BF6 - 0.037 * BF8 ...
- 0.241 * BF10 + 0.018 * BF11 - 0.008 * BF12 ...
+ 0.059 * BF13 - 0.144 * BF17 - 0.116 * BF18 ...
+ 0.010 * BF19 - 0.012 * BF20 + 0.085 * BF22 ...
+ 0.011 * BF23 + 0.049 * BF25 - 0.159 * BF26 ...
- 0.016 * BF28 - 0.138 * BF29 + 0.010 * BF31 ...
+ 0.016 * BF35 + 0.018 * BF36 + 0.246 * BF37 ...
- 0.417 * BF38 + 0.052 * BF39 - 0.005 * BF41 ...
+ 0.021 * BF43 + 0.006 * BF44 - 0.047 * BF45 ...
- 0.051 * BF47 - 0.014 * BF48 - 0.113 * BF50 ...
+ 0.019 * BF51 + 0.007 * BF53 + 0.017 * BF54 ...
- 0.007 * BF56 - 0.098 * BF58 + 0.011 * BF59 ...
- 0.016 * BF62 - 0.012 * BF63 + 0.113 * BF65 ...
+ 0.016 * BF66 + 0.040 * BF67 - 0.065 * BF69 ...
- 0.018 * BF71 + 0.014 * BF73 - 0.009 * BF75 ...
- 0.008 * BF76 - 0.032 * BF78 + 0.032 * BF79 ...
+ 0.011 * BF80;
function Y= ovl_mars( FWSEG_VA, V5, V6, V7, V8, V9, V10, V11, V12, ...
V13, V14, V15, V16, V17, V18, V19, V20, ...
V21, V22, V23, V24, V25, V26, V27, V28)
BF1 = max(0, V21 - 4.671);
BF3 = max(0, V6 - 5.396);
BF4 = max(0, 5.396 - V6 );
BF7 = max(0, V11 - 7.884);
BF8 = max(0, 7.884 - V11 );
BF9 = max(0, FWSEG_VA + 7.229) * BF1;
BF10 = max(0, - 7.229 - FWSEG_VA ) * BF1;
BF11 = max(0, V19 - 8.128) * BF1;
BF12 = max(0, 8.128 - V19 ) * BF1;
BF13 = max(0, V28 - 7.918);
BF14 = max(0, 7.918 - V28 );
BF15 = max(0, V5 + 2.888) * BF14;
BF16 = max(0, - 2.888 - V5 ) * BF14;
BF17 = max(0, V24 - 2.924) * BF8;
BF18 = max(0, 2.924 - V24 ) * BF8;
BF20 = max(0, 9.071 - V16 ) * BF15;
BF21 = max(0, V10 - 6.286) * BF14;
BF22 = max(0, 6.286 - V10 ) * BF14;
BF24 = max(0, V23 - 5.173);
BF25 = max(0, 5.173 - V23 );
BF26 = max(0, V26 - 8.987);
BF29 = max(0, 12.216 - V27 ) * BF3;
BF30 = max(0, V8 - 4.306) * BF16;
BF34 = max(0, V23 - 7.630) * BF21;
BF35 = max(0, 7.630 - V23 ) * BF21;
BF37 = max(0, 3.638 - V7 ) * BF1;
BF39 = max(0, 8.337 - V21 ) * BF17;
BF41 = max(0, 1.590 - V5 ) * BF11;
BF43 = max(0, 13.993 - V8 ) * BF11;
BF44 = max(0, V14 - 5.993) * BF25;
BF45 = max(0, 5.993 - V14 ) * BF25;
BF46 = max(0, V24 - 1.035);
BF47 = max(0, 1.035 - V24 );
BF49 = max(0, 8.915 - V23 ) * BF12;
BF51 = max(0, - 0.004 - FWSEG_VA );
BF52 = max(0, V27 - 6.520) * BF24;
BF53 = max(0, 6.520 - V27 ) * BF24;
BF54 = max(0, V7 - 11.484) * BF8;
BF55 = max(0, 11.484 - V7 ) * BF8;
BF57 = max(0, 5.742 - V17 ) * BF25;
BF58 = max(0, V12 - 6.949) * BF12;
BF59 = max(0, 6.949 - V12 ) * BF12;
BF60 = max(0, V25 - 9.203) * BF45;
BF63 = max(0, 1.887 - V13 ) * BF7;
BF65 = max(0, 9.498 - V26 ) * BF15;
BF66 = max(0, V5 - 6.566) * BF22;
BF71 = max(0, 13.239 - V19 ) * BF46;
BF72 = max(0, V19 - 9.925) * BF55;
BF77 = max(0, 3.430 - V22 ) * BF18;
BF78 = max(0, V27 - 6.513) * BF45;
BF79 = max(0, 6.513 - V27 ) * BF45;
BF81 = max(0, 12.511 - V18 );
BF82 = max(0, V11 - 6.777) * BF81;
BF83 = max(0, 6.777 - V11 ) * BF81;
BF85 = max(0, 3.433 - V5 ) * BF47;
BF87 = max(0, - 3.524 - FWSEG_VA ) * BF47;
BF88 = max(0, V27 - 11.604) * BF9;
BF91 = max(0, 8.845 - V26 ) * BF52;
BF92 = max(0, V14 - 5.931) * BF82;
BF93 = max(0, 5.931 - V14 ) * BF82;
BF94 = max(0, V21 - 7.245) * BF25;
BF95 = max(0, 7.245 - V21 ) * BF25;
BF96 = max(0, V14 - 5.323) * BF7;
BF98 = max(0, V10 - 6.248) * BF71;
BF100 = max(0, V18 - 0.602) * BF95;
Y = 2.936 + 0.047 * BF1 + 0.061 * BF3 - 0.084 * BF4 - 0.139 * BF8 ...
- 0.064 * BF10 - 0.030 * BF12 - 0.103 * BF13 ...
- 0.039 * BF14 + 0.020 * BF17 - 0.002 * BF20 ...
- 0.005 * BF22 - 0.114 * BF25 - 0.090 * BF26 ...
- 0.011 * BF29 + 0.010 * BF30 + 0.009 * BF34 ...
+ 0.002 * BF35 + 0.079 * BF37 - 0.006 * BF39 ...
+ 0.007 * BF41 - 0.003 * BF43 + 0.017 * BF44 ...
+ 0.076 * BF47 + 0.009 * BF49 + 0.016 * BF51 ...
- 0.042 * BF53 - 0.079 * BF54 - 0.030 * BF57 ...
- 0.018 * BF58 - 0.009 * BF59 - 0.119 * BF60 ...
- 0.210 * BF63 - .456802E-03 * BF65 + 0.028 * BF66 ...
+ 0.020 * BF72 + 0.011 * BF77 + 0.005 * BF78 ...
+ 0.003 * BF79 - 0.049 * BF81 + 0.012 * BF83 ...
- 0.030 * BF85 + 0.070 * BF87 + 0.008 * BF88 ...
- 0.008 * BF91 + 0.010 * BF92 + 0.003 * BF93 ...
+ 0.022 * BF94 - 0.038 * BF96 + .933766E-03 * BF98 ...
+ 0.002 * BF100;
function distortion = fwseg(clean_speech, processed_speech,sample_rate)
% ----------------------------------------------------------------------
% Check the length of the clean and processed speech. Must be the same.
% ----------------------------------------------------------------------
clean_length = length(clean_speech);
processed_length = length(processed_speech);
if (clean_length ~= processed_length)
disp('Error: Files must have same length.');
return
end
% ----------------------------------------------------------------------
% Global Variables
% ----------------------------------------------------------------------
winlength = round(30*sample_rate/1000); % window length in samples
skiprate = floor(winlength/4); % window skip in samples
max_freq = sample_rate/2; % maximum bandwidth
num_crit = 25; % number of critical bands
n_fft = 2^nextpow2(2*winlength);
n_fftby2 = n_fft/2; % FFT size/2
% ----------------------------------------------------------------------
% Critical Band Filter Definitions (Center Frequency and Bandwidths in Hz)
% ----------------------------------------------------------------------
cent_freq(1) = 50.0000; bandwidth(1) = 70.0000;
cent_freq(2) = 120.000; bandwidth(2) = 70.0000;
cent_freq(3) = 190.000; bandwidth(3) = 70.0000;
cent_freq(4) = 260.000; bandwidth(4) = 70.0000;
cent_freq(5) = 330.000; bandwidth(5) = 70.0000;
cent_freq(6) = 400.000; bandwidth(6) = 70.0000;
cent_freq(7) = 470.000; bandwidth(7) = 70.0000;
cent_freq(8) = 540.000; bandwidth(8) = 77.3724;
cent_freq(9) = 617.372; bandwidth(9) = 86.0056;
cent_freq(10) = 703.378; bandwidth(10) = 95.3398;
cent_freq(11) = 798.717; bandwidth(11) = 105.411;
cent_freq(12) = 904.128; bandwidth(12) = 116.256;
cent_freq(13) = 1020.38; bandwidth(13) = 127.914;
cent_freq(14) = 1148.30; bandwidth(14) = 140.423;
cent_freq(15) = 1288.72; bandwidth(15) = 153.823;
cent_freq(16) = 1442.54; bandwidth(16) = 168.154;
cent_freq(17) = 1610.70; bandwidth(17) = 183.457;
cent_freq(18) = 1794.16; bandwidth(18) = 199.776;
cent_freq(19) = 1993.93; bandwidth(19) = 217.153;
cent_freq(20) = 2211.08; bandwidth(20) = 235.631;
cent_freq(21) = 2446.71; bandwidth(21) = 255.255;
cent_freq(22) = 2701.97; bandwidth(22) = 276.072;
cent_freq(23) = 2978.04; bandwidth(23) = 298.126;
cent_freq(24) = 3276.17; bandwidth(24) = 321.465;
cent_freq(25) = 3597.63; bandwidth(25) = 346.136;
bw_min = bandwidth (1); % minimum critical bandwidth
% ----------------------------------------------------------------------
% Set up the critical band filters. Note here that Gaussianly shaped
% filters are used. Also, the sum of the filter weights are equivalent
% for each critical band filter. Filter less than -30 dB and set to
% zero.
% ----------------------------------------------------------------------
min_factor = exp (-30.0 / (2.0 * 2.303)); % -30 dB point of filter
for i = 1:num_crit
f0 = (cent_freq (i) / max_freq) * (n_fftby2);
all_f0(i) = floor(f0);
bw = (bandwidth (i) / max_freq) * (n_fftby2);
norm_factor = log(bw_min) - log(bandwidth(i));
j = 0:1:n_fftby2-1;
crit_filter(i,:) = exp (-11 *(((j - floor(f0)) ./bw).^2) + norm_factor);
crit_filter(i,:) = crit_filter(i,:).*(crit_filter(i,:) > min_factor);
end
% ----------------------------------------------------------------------
% For each frame of input speech, calculate the Weighted Spectral
% Slope Measure
% ----------------------------------------------------------------------
num_frames = floor(clean_length/skiprate-(winlength/skiprate)); % number of frames
start = 1; % starting sample
window = 0.5*(1 - cos(2*pi*(1:winlength)'/(winlength+1)));
distortion=zeros(num_frames,num_crit);
for frame_count = 1:num_frames
% ----------------------------------------------------------
% (1) Get the Frames for the test and reference speech.
% Multiply by Hanning Window.
% ----------------------------------------------------------
clean_frame = clean_speech(start:start+winlength-1);
processed_frame = processed_speech(start:start+winlength-1);
clean_frame = clean_frame.*window;
processed_frame = processed_frame.*window;
% ----------------------------------------------------------
% (2) Compute the magnitude Spectrum of Clean and Processed
% ----------------------------------------------------------
clean_spec = abs(fft(clean_frame,n_fft));
processed_spec = abs(fft(processed_frame,n_fft));
% normalize so that spectra have unit area ----
clean_spec=clean_spec/sum(clean_spec(1:n_fftby2));
processed_spec=processed_spec/sum(processed_spec(1:n_fftby2));
% ----------------------------------------------------------
% (3) Compute Filterbank Output Energies
% ----------------------------------------------------------
clean_energy=zeros(1,num_crit);
processed_energy=zeros(1,num_crit);
error_energy=zeros(1,num_crit);
for i = 1:num_crit
clean_energy(i) = sum(clean_spec(1:n_fftby2) ...
.*crit_filter(i,:)');
processed_energy(i) = sum(processed_spec(1:n_fftby2) ...
.*crit_filter(i,:)');
error_energy(i)=max((clean_energy(i)-processed_energy(i))^2,eps);
end
SNRlog=10*log10((clean_energy.^2)./error_energy);
distortion(frame_count,:)=min(max(SNRlog,-10),35);
start = start + skiprate;
end