%lossyFLAC.m % David Robinson, 2007 % This code is open source. I'll pick a licence later (need advice), but: % Any changes to this MATLAB code must be shared freely % Any changes to the algorithm must be shared freely % Any other implementations / interpretations must be shared freely % You are free to include any implementation in commercial and/or closed source code, but must give an acknowledgement in help/about or similar. % No warrantee / guarantee. This is work in progress and not debugged. %=== % external variables (i.e. 'switches' when implemented properly!) foldername='D:\audio\lossless test\problem samples\'; filename='Image32k'; % don't put .wav at end noise_threshold_shift=0; % noise threshold shift. average level of added quantisation % noise relative to the lowest amplitude frequency bin (default=0=equal!) low_frequency_limit=20; high_frequency_limit=16000; % Frequency range over which to check for lowest amplitude signal fix_clipped=0; %=== % internal constants (i.e. no need to change them) flac_blocksize=4096; noise_averages=1000; inaudible=1E-10; % load file [inaudio,fs,bs]=waveread([foldername filename '.wav']); % This pads it to 23 bits: It's treated as 24 bit data at the end, effectively dropping it by 6dB if fix_clipped==1, bs=23; end [samples channels]=size(inaudio); % create integer copy (MATLAB wav(e)read loads audio with the range +/-1) inaudio_int=inaudio.*(2^(bs-1))+inaudible; % Set up the FFT analysis lengths - define these in time (seconds) % (will then use the nearest power of two based on sampling frequency) clear analysis_time spreading_function fft_length low_frequency_bin high_frequency_bin reference_thresholds reference_threshold min_bin max_bin bits_to_remove bits_to_remove_table analysis_time(1)=2.0E-2; % 20ms analysis_time(2)=1.5E-3; % 1.5ms %analysis_time(2)=3.0E-3; number_of_analyses=length(analysis_time); spreading_function{1}=[0.25,0.25,0.25,0.25]; spreading_function{2}=[0.25,0.25,0.25,0.25]; % spreading function to apply to FFT before determining lowest amplitude % Keep peak at centre, even if it means padding with zeros % Loop through each analysis length (typically only two) % and set the FFT values for analysis_number=1:number_of_analyses, % Calculate the closest FFT length fft_length(analysis_number)=min(2^round(log10(analysis_time(analysis_number)*fs)/log10(2)),flac_blocksize); % Calculate which FFT bin corresponds to the low frequency limit low_frequency_bin(analysis_number)=round(fft_length(analysis_number)*low_frequency_limit/fs+((length(spreading_function)-1)/2)); if low_frequency_bin(analysis_number)<2, low_frequency_bin(analysis_number)=2; end if low_frequency_bin(analysis_number)>fft_length(analysis_number)/2, error('low frequency too high'); end % Calculate which FFT bin corresponds to the high frequency limit high_frequency_bin(analysis_number)=round(fft_length(analysis_number)*high_frequency_limit/fs+((length(spreading_function)-1)/2)); if high_frequency_bin(analysis_number)<2, error('high frequency too low'); end if high_frequency_bin(analysis_number)>fft_length(analysis_number)/2, high_frequency_bin(analysis_number)=fft_length(analysis_number)/2; end end variables_filename=['lossy_variables__fs' num2str(fs) '_bs' num2str(bs) '_noa' num2str(number_of_analyses) '_fft' num2str(fft_length) '_lfb' num2str(low_frequency_bin) '_hfb' num2str(high_frequency_bin) '_nts' num2str(noise_threshold_shift) '.m']; % Find out if we've stored the quantisation noise thresholds before if exist(variables_filename,'file'), load('-MAT',variables_filename) % If not, calculate quantisation noise at each bit in these FFTs else for analysis_number=1:number_of_analyses, reference_thresholds(1:noise_averages,1:bs)=0; for av_number=1:noise_averages, noise_sample=rand(fft_length(analysis_number),1); for bits_to_remove=1:bs, % This models the quantisation noise introduced by truncating % the last "bits_to_remove" bits from the audio data this_noise_sample=floor(noise_sample.*((2^bits_to_remove)))-(2^(bits_to_remove-1)); fft_result=20*log10(conv(abs(fft(this_noise_sample.*hanning(fft_length(analysis_number)))),spreading_function{analysis_number})); reference_thresholds(av_number,bits_to_remove)=mean(fft_result(low_frequency_bin(analysis_number):high_frequency_bin(analysis_number))); end end reference_threshold{analysis_number}=mean(reference_thresholds)-noise_threshold_shift; for threshold=1:round(20*log10(2^(bs+4))), if isempty(find(reference_threshold{analysis_number}samples, block_end=samples; end % Don't jump past end of file! for analysis_number=1:number_of_analyses, first_block=(block_start-1)/(fft_length(analysis_number)/2); last_block=first_block+(flac_blocksize/(fft_length(analysis_number)/2)); if first_block<1, first_block=1; end % Don't jump before start of file if last_block>min_bin_length(analysis_number), last_block=min_bin_length(analysis_number); end % Don't jump past end of file! for channel=1:channels, this_min_bin=round(min(min_bin{analysis_number}(first_block:last_block,channel))); if this_min_bin<1, % i.e. if it's quieter than quantisation noise at the least significant bit bits_to_remove_table(analysis_number,channel)=0; % don't remove any bits! else bits_to_remove_table(analysis_number,channel)=threshold_index{analysis_number}(this_min_bin); end end end bits_to_remove(block_number)=min(min(bits_to_remove_table)); if bits_to_remove(block_number)>0, inaudio_int(block_start:block_end,1:channels)=round(inaudio_int(block_start:block_end,1:channels)/(2^bits_to_remove(block_number))).*2^bits_to_remove(block_number); end end if fix_clipped==1, bs=24; end outaudio=inaudio_int./(2^(bs-1)); wavewrite(outaudio,fs,bs,[foldername filename ' lossy.wav'])