www.pudn.com > MATLAB.rar > noise_inv_f.m, change:2010-01-07,size:2066b

function hout=noise_inv_f(L); 
% L= filter length/2 
% Create a filter kernel (FIR coefficients) that approximates a 1/f frequency response. 
% Dick Benson  September 2004 
% Copywrite 2004-2010 The MathWorks, Inc. 
if nargin==0 
   L     = 128;             % filter kernel length/2 
f     =(0:(L-1))/L;     % frequency vector normalized to Fs/2 
shape = (1./(f+0.1/L)).^0.5;  % desired shape of spectrum (note +0.1/L to limit infinity ) ,  
                              % added .^0.5 to make POWER drop as 1/f 
shape = shape/shape(1); % normalize DC gain to 1 (0dB) 
ph    = pi*(0:L-1);     % phase shift to center impulse response in time window 
% Need to create a complex spectral description  
% which, when transformed with the ifft, 
% creates a REAL impulse response. 
mag=[shape,0, fliplr(shape(2:end))]; % construct magnitude 
phase=[-ph,0,fliplr(ph(2:end))];     % and phase  
H=mag.*exp(j*phase);  % complex representation 
h=ifft(H);            % impulse reponse  
%  sanity checks ....  
plot((h)); title('1/f filter impulse response (no window)'); 
% If the filter is used as is, there will be ripples in the frequency reponse,  
% due to the abrupt transistions at the impulse response endpoints. 
% To mitigate this, apply a Hanning window and then compare final response with 
% desired.  
H2=fft(hout);            % take fft of windowed filter kernel 
% hf=figure; 
title('overlayed freq responses of desired and attained 1/f shape'); 
xlabel('Normalized to Fs'); ylabel('dB') 
% Finally, scale the filter so that 1 Watt of (broad band) input noise  
% produces 1 Watt of output noise.  
df   = f(2)-f(1); 
Pout_T = sum(hout.^2);       % power out for 1 W input 
hout = hout/sqrt(Pout_T);    % use filter_check.mdl to confirm power level and shape 
f_inv = hout; 
save f_inv f_inv;            % save for use by other examples.