Why does pwelch signal magnitude change with frequency?

I am sampling a signal with high and low frequency components that I know should occur. My sampling rate is 2.56*10^4 Hz, and I have 577664 data points in the input vector x.

I am using pwelch and a blackman-harris window with 50% overlap to estimate the magnitude of my signal in the frequency domain with the code example below.

When I use the code below, the low frequency magnitudes (0 to 200 Hz) are about 30 dB greater than what I expect and the high frequency components (600 to 1000 Hz) are about 20 dB lower than I expect. Why would this happen and how can I fix it?

Thank you.

Code Example:

    fs = 2.56*10^4;
    n = fs/2;

    win = blackmanharris(floor(n)); % Window Function

    [magFFTX,f] = pwelch(x,win,floor(n/2),[],fs); % Calculate Mag and f

    fbin = f(2,1) - f(1,1); % Calculate frequency bin

    %CG = 0.5400; % Coherent Gain of Hamming Window
    %NG = 0.3974; % Noise Gain of Hamming Window

    %CG = 0.5000; % Coherent Gain of Hanning Window
    %NG = 0.3750; % Noise Gain of Hanning Window

    CG = 0.3587; % Coherent Gain of Blackman-Harris Window
    NG = 0.2580; % Noise Gain of Blackman-Harris Window

    scalingFactor = NG.*fbin./(CG.^2);

    magFFTX = magFFTX.*scalingFactor;

    magFFTX = magFFTX ./(10^-6); % Convert to micro g
    magFFTX = 20.*log10(magFFTX ); % Convert to dB

    plot(f,magFFTX);
    xlim([0 1000]);

[Image Analyst]

You forgot to supply x. I'm not doing anything until I get x.

[John BG]

IA, hugh didn't forget the signal: it's a filter therefore it should work with any signal you plug in.

Have a go with following:

t=[0:1/fs:1-1/fs];x=10*sin(2*pi*1e3*t)

try different frequencies to reproduce what is mentioned in the question

[Image Analyst]

Well he did forget the signal, x, because when you run his code, simply copying and pasting, it says this:

Undefined function or variable 'x'.
Error in test3 (line 6)
    [magFFTX,f] = pwelch(x,win,floor(n/2),[],fs); % Calculate Mag and f

John took a stab at creating and supplying an x and got a signal of 25,600 elements long. His actual signal is 577,664 elements long and is missing. Wouldn't it be so much better for posters to supply their data than requiring us to invent/create/synthesize data just to be able to run their code?

[J H]

Here is a signal that replicates what I am talking about:

The magnitude is wrong in both linear and log space.

Could it be due to high frequency noise getting aliased? If so, how can I filter that out properly with a lowpass filter?

Thanks.

t = (0:1/fs:1-1/fs);
x = 2*sin(2*pi*97*t) + 1.5*sin(2*pi*149*t) + 1*sin(2*pi*213*t) + 
0.75*sin(2*pi*416*t) + 0.5*sin(2*pi*555*t)+ 0.25*sin(2*pi*611*t)+ 0.1*sin(2*pi*675*t)+ 0.1*sin(2*pi*700*t)+ 0.1*sin(2*pi*760*t)+ 0.1*sin(2*pi*813*t);
x = x + 10.*randn(size(t)); % Add Some Noise

[John BG]

the amplitude of the MATLAB fft function varies with the amount of time samples. For exactly same periodic signal, the y of the question:

y1=y([1:12]);figure(2);plot(abs(fft(y1))); hold all
y2=y([1:24]);figure(2);plot(abs(fft(y2)))
y3=y([1:48]);figure(2);plot(abs(fft(y3)))

and IA, no one is, has or will take stabs here, this is just an online forum: people write things, right or wrong, just words, just software no hardware.

[J H]

Thanks.

So in general, what could cause lower frequencies to have higher amplitudes and higher frequencies to have lower amplitudes (than my reference)?

I do not think that it is aliasing because that would raise the amplitude on all frequencies. I also tried to use MATLABs lowpass filter to filter out everything above the nyquist frequency and that did not change my result.

I doubt it is spectral leakage because that would lower the amplitude on all frequencies by spreading out the energy.

Am I forgetting to normalize my output by some scaling factor that varies with frequency?

Thanks again!

[John BG]

it's the Barlett window used in pwelch deforming the measurement over a wide frequency range, let me explain:

1. from MATLAB help

Welch's Overlapped Segment Averaging Spectral Estimation

The periodogram is not a consistent estimator of the true power spectral density of a wide-sense stationary process.

Welch's technique to reduce the variance of the periodogram breaks the time series into segments, usually overlapping.

Welch's method computes a modified periodogram for each segment and then averages these estimates to produce the estimate of the power spectral density. Because the process is wide-sense stationary and Welch's method uses PSD estimates of different segments of the time series, the modified periodograms represent approximately uncorrelated estimates of the true PSD and averaging reduces the variability.

2. from

http://www.dcs.warwick.ac.uk/~feng/teaching/PowerSpectrum.pdf

point 12

3. With Uniform, no windowing, there is no such spectral deformation.

if you find these lines useful would you please mark my answer as Accepted Answer?

To any other reader, if you find this answer of any help please click on the thumbs-up vote link,

thanks in advance for time and attention

John BG

<mailto:[e-mail] [e-mail]>

[Image Analyst]

John - he can't mark it as accepted since you didn't leave an answer - you left a comment. If you move your good Comments into an Answer below, I'd vote for it and he'd also be able to Accept it and Vote for it.

[Alexander Tsarev]

.

[J H]

Thanks John.

One last question.

I have to use windowing because my signal is not completely uniform (like a code generate sine wave). It does have some noise and nonuniformity to it because its real sampled data.

So, how do signal analysts normally use windows to deal with discontinuities, while preserving the true power spectral density i.e. is there some different Matlab function or Fourier definition that I can use in my case?

Thanks again!

[John BG]

ok,

1.

is the comb of tones you supplied as example, just an example, or is it the real signal you want to analyse?

I ask this because the window to apply depends upon the type of signal you want to analyse.

2.

would you consider splitting the input spectrum into smaller bands and analyse them separately?

3.

If the supplied signal is not the signal you want to analyse, would it be possible to have a look at the the actual signal?

4.

what kind of signals are you working with

Are you analysing coupled interference across different loops?

or are you monitoring wireless traffic picked by a mobile communications tower?

[J H]

1. The comb of tones I supplied as an example is just an example.

2. I would consider this, if it may help. Would I split the data before FFT, and how/where should I split it?

3. I am not allowed to provide the actual signal. Sorry.

4. The signal comes from accelerometers, placed on a vibrating cable. The cable is forced by a fluid flowing at a certain velocity past it. I am taking the maximum at each timestep of specifically placed accelerometer data to get the natural frequencies of the cable. This is the power spectrum that I am trying to plot.

Thanks.