DSPRelated.com
Code

Frequency Response Creator

Ron April 11, 2011 Coded in Matlab

To be used with a sweeping tone, this program takes the output and calculates the magnitude response at each frequency.

This is useful for determining the frequency response of any audio device that you are testing. Be sure to set the sampling frequency.

[header, s] = hdrload('data_file.txt');

num = 2;
x = zeros(ceil(44100/num),1);                 %channel 1
y = zeros(ceil(44100/num),1);                 %channel 2
 
Fs = 44100;                     % Sampling frequency
T = 1/Fs;                       % Sample time
L = Fs/num;                  	  % Length of window
g = floor(num*length(s)/Fs);    % number of windows
H = 0;
count = 50;
max_index = 1;                  % pointer for max value
max = 0;                        % max value
max_index2 = 1;                 % pointer for max value
max2 = 0;                       % max value
 
z = zeros(g,2);
z2 = zeros(g,2);
 
for i = 1:1:g-3
    
% creating windows for channel one and two
    for j = 1:1:L
       k = j + L*i;
       x(j,1) = s(k+i,1);
       y(j,1) = s(k+i,2);
    end
 
NFFT = 2^nextpow2(L);            % Next power of 2 from length of y
YY = fft(y,NFFT)/L;
Y = fft(x,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1); % creating frequency vector
Y2 = 2*abs(Y(1:NFFT/2+1));
YY2 = 2*abs(YY(1:NFFT/2+1));
 
%%channel 1
max = Y(max_index,1);
 
%%channel 2 
if(max_index2 < length(YY2))
max2 = YY(max_index2,1);
 
b = length(YY2) - max_index2 - 2*count;
if b > 0
    a = max_index2 + 2*count; 
else
    a = length(YY2);
end
 
for j = max_index2:1:(a-1)
    if max2 < YY2(j+1,1)
        max2 = YY2(j+1,1);
        max = Y2(j+1,1);
        max_index2 = j+1;
    end
end
end
    z2(i+1,1) = abs(max2);
    z2(i+1,2) = f(1,max_index2);
    z(i+1,1) = abs(max);
    z(i+1,2) = f(1,max_index2);
    if(max_index2 > 4*count)
    max_index2 = max_index2 - count;
    end

% RECORD MOVIE
% semilogy(z2(:,2),z2(:,1),f,YY2,z(:,2),z(:,1),f,Y2);
% axis([0 20000 10^-3 Inf]);

end
 
semilogy(z2(:,2),z2(:,1),z(:,2),z(:,1));
axis([0 20000 0.00001 Inf]);