function [ f, p, pint, pintfx ] = avg_spectrum(y, fx) % AVG_SPECTRUM average power spectrum. Uses a Hodie window. % Inputs: % Y K x N matrix. Each of K columns is interpreted as a separate signal. % A hodie window is applied to each column and the K spectra are averaged. % FX Frequency of fundamental, relative to sampling rate (in range 0 ... % 0.5). Used to calculate PINTFX. Optional parameter. % Outputs: % F Frequency vector, 0 ... 0.5, length N/2 + 1. % P Averaged power spectrum, length N/2 [dBWN]. % PINT Integrated power spectrum, length N/2. PINT(end) is the power of the % signal. % PINTFX Same as PINT, except that the signal power is excluded from the % integral. Computed only when fx is specified. [ N, K ] = size(y); % frequency vector f = linspace(0, 0.5, N/2+1)'; % window w = hodiewindow(N); w = w*ones(1, K); % average magnitude from 0 to fs/2 p = fft(y .* w); % windowed spectrum p(end/2+2:end, :) = []; % frequencies 0 to fs/2 p = p .* conj(p); % power if K > 1 p = sum(p')' / K; % average end p(1) = 0.5*p(1); % dc fix p(end) = 0.5*p(end); p = p/N; % noise power integral pint = 2/N*cumsum(p); % noise integral if nargin > 1 nx = round(N * fx + 1); % fundamental center bin pp = 2/N*p; x1 = max([1, nx-7]); x2 = min([length(pp), nx+7]); pp(x1:x2) = 0; % zero out fundamental pp(1) = 0; % zero out dc term pintfx = cumsum(pp); end