image thumbnail

Noise Measuremet with PCSGU250 and Matlab Implementation

by

 

16 Jul 2012 (Updated )

Noise measurement using data from the Velleman PCSGU250 measurement complex.

filterA(x, fs, plotFilter)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         Generation of an A-weighting filter          %
%              with Matlab Implementation              %
%                                                      %
% Authors: Douglas R. Lanman, 11/21/05                 %
%          Hristo Zhivomirov, 01/07/14                 %
%                                                      %
% Define filter coefficients:                          %
% IEC 61672:2003                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function xA = filterA(x, fs, plotFilter)

% function: xA = filterA(x, fs, plotFilter)
% x - singnal in the time domain
% fs - sampling frequency, Hz
% type 1 in the place of plotFilter if one want a plot of freq response
% xA - filtered singnal in the time domain (row-vector)

% filter coefficients
c1 = 3.5041384e16;
c2 = 20.598997^2;
c3 = 107.65265^2;
c4 = 737.86223^2;
c5 = 12194.217^2;

% represent x as row-vector if it is not
if size(x,1) > 1
    x = x';
end

% calculate xlen
xlen = length(x);

% calculate the number of unique points
NumUniquePts = ceil((xlen+1)/2);

% take fft of x
X = fft(x);

% fft is symmetric, throw away second half
X = X(1:NumUniquePts);

% frequency vector with NumUniquePts points
f = (0:NumUniquePts-1)*fs/xlen;

% evaluate A-weighting filter
f = f.^2;
num = c1*f.^4;
den = ((c2+f).^2) .* (c3+f) .* (c4+f) .* ((c5+f).^2);
A = num./den;

% filtering
XA = X.*A;

% perform ifft
if rem(xlen, 2)                     % odd xlen excludes Nyquist point
    % reconstruct the whole spectrum
    XA = [XA conj(XA(end:-1:2))];
    
    % take ifft of XA
    xA = real(ifft(XA));
else                                % even xlen includes Nyquist point
    % reconstruct the whole spectrum
    XA = [XA conj(XA(end-1:-1:2))];
    
    % take ifft of XA
    xA = real(ifft(XA));
end

% plot A-weighting filter (if enabled)
if exist('plotFilter') && plotFilter
    
    % plot using dB scale
    figure
    semilogx(sqrt(f), 10*log10(A))
    grid on
    set(gca, 'FontName', 'Times New Roman', 'FontSize', 13)
    title('A-weighting Filter')
    xlabel('Frequency, Hz')
    ylabel('Magnitude, dB')
        
    % plot using linear scale
    figure
    plot(sqrt(f)/1000, A)
    grid on
    set(gca, 'FontName', 'Times New Roman', 'FontSize', 13)
    title('A-weighting Filter')
    xlabel('Frequency, kHz')
    ylabel('Amplitude')
        
end
end

Contact us