Thread Subject:

PSD and windowing function

I have a query concerning PSD. I am trying to use MATLAB for this, but the dominant frequencies I am getting are very scattered and I need to incorporate a windowing function to smooth out all the peaks. I would like to incorporate to begin with a simple cosine windowing function of the form w(kdt)=cos (Pi*kdt/2 *t). where dt is resolution of the time delay. I am unsure how to put this into my program, and this is simpler than the window function inbuilt in Matlab. The program is as follows (H1 in this example has 30000 data points):

H1=a(:,2);
t=0.001:0.001:30;
time_interval=t(2)-t(1);
sizeHl=30000;
Y = fft(Hl);
fr = (1/time_interval)*(1:15000)/sizeHl;
Pyy = Y.* conj(Y)/sizeHl;
Pyy_interval=Pyy(2:15001);
newcounter=1;
 
[C,I] = max(Pyy_interval);
f= (I/time_interval)/sizeHl;

frequency =f'

Any suggestions to how to improve the above program would be much appreciated.
Many thanks, Safa

On Aug 30, 9:50 am, "Safa " <enxs...@nottingham.ac.uk> wrote:
> I have a query concerning PSD. I am trying to use MATLAB for this, but the dominant frequencies I am getting are very scattered and I need to incorporate a windowing function to smooth out all the peaks. I would like to incorporate to begin with a simple cosine windowing function of the form w(kdt)=cos (Pi*kdt/2 *t). where dt is resolution of the time delay. I am unsure how to put this into my program, and this is simpler than the window function inbuilt in Matlab. The program is as follows (H1 in this example has 30000 data points):
>
> H1=a(:,2);
> t=0.001:0.001:30;
> time_interval=t(2)-t(1);
> sizeHl=30000;
> Y = fft(Hl);
> fr = (1/time_interval)*(1:15000)/sizeHl;
> Pyy = Y.* conj(Y)/sizeHl;
> Pyy_interval=Pyy(2:15001);
> newcounter=1;
>
> [C,I] = max(Pyy_interval);
> f= (I/time_interval)/sizeHl;
>
> frequency =f'
>
> Any suggestions to how to improve the above program would be much appreciated.
> Many thanks, Safa

Windowing the data will not make the PSD smoother, it's purpose is to
remove the end effects from the time series.
To smooth the spectrum, you need to either do ensemble averaging
(split the record into blocks, calculate the PSD of each and average),
or frequency averaging (accumulate the energy into large frequency
blocks), or both.

BTW, your Pyy_interval is not a PSD - it must be divided by the
frequency interval, ie Pyy_interval*time_interval*sizeHI.

Also, how come H1 has changed to HI after line 1 in your code?

On Aug 30, 1:04 am, TideMan <mul...@gmail.com> wrote:
> On Aug 30, 9:50 am, "Safa " <enxs...@nottingham.ac.uk> wrote:
>
>
>
>
>
> > I have a query concerning PSD. I am trying to use MATLAB for this, but the dominant frequencies I am getting are very scattered and I need to incorporate a windowing function to smooth out all the peaks. I would like to incorporate to begin with a simple cosine windowing function of the form w(kdt)=cos (Pi*kdt/2 *t). where dt is resolution of the time delay. I am unsure how to put this into my program, and this is simpler than the window function inbuilt in Matlab. The program is as follows (H1 in this example has 30000 data points):
>
> > H1=a(:,2);
> > t=0.001:0.001:30;
> > time_interval=t(2)-t(1);
> > sizeHl=30000;
> > Y = fft(Hl);
> > fr = (1/time_interval)*(1:15000)/sizeHl;
> > Pyy = Y.* conj(Y)/sizeHl;
> > Pyy_interval=Pyy(2:15001);
> > newcounter=1;
>
> > [C,I] = max(Pyy_interval);
> > f= (I/time_interval)/sizeHl;
>
> > frequency =f'
>
> > Any suggestions to how to improve the above program would be much appreciated.
> > Many thanks, Safa
>
> Windowing the data will not make the PSD smoother,

It will.

> it's purpose is to
> remove the end effects from the time series.

Its purpose is to reduce the variance of the periodogram.

> To smooth the spectrum, you need to either do ensemble averaging
> (split the record into blocks, calculate the PSD of each and average),
> or frequency averaging (accumulate the energy into large frequency
> blocks), or both.

If you do the maths, you might find that windowing in time
domain is equivalent to convolution in frequency domain.
Convolution is, in turn, equivalent to weighted averaging
of the sequence, in this case the spectrum.

Rune

Rune Allnor <allnor@tele.ntnu.no> wrote in message <b4b62db2-03af-4a5d-9e48-5317a5a72b2f@n3g2000yqb.googlegroups.com>...
> On Aug 30, 1:04 am, TideMan <mul...@gmail.com> wrote:
> > On Aug 30, 9:50 am, "Safa " <enxs...@nottingham.ac.uk> wrote:
> >
> >
> >
> >
> >
> > > I have a query concerning PSD. I am trying to use MATLAB for this, but the dominant frequencies I am getting are very scattered and I need to incorporate a windowing function to smooth out all the peaks. I would like to incorporate to begin with a simple cosine windowing function of the form w(kdt)=cos (Pi*kdt/2 *t). where dt is resolution of the time delay. I am unsure how to put this into my program, and this is simpler than the window function inbuilt in Matlab. The program is as follows (H1 in this example has 30000 data points):
> >
> > > H1=a(:,2);
> > > t=0.001:0.001:30;
> > > time_interval=t(2)-t(1);
> > > sizeHl=30000;
> > > Y = fft(Hl);
> > > fr = (1/time_interval)*(1:15000)/sizeHl;
> > > Pyy = Y.* conj(Y)/sizeHl;
> > > Pyy_interval=Pyy(2:15001);
> > > newcounter=1;
> >
> > > [C,I] = max(Pyy_interval);
> > > f= (I/time_interval)/sizeHl;
> >
> > > frequency =f'
> >
> > > Any suggestions to how to improve the above program would be much appreciated.
> > > Many thanks, Safa
> >
> > Windowing the data will not make the PSD smoother,
>
> It will.
>
> > it's purpose is to
> > remove the end effects from the time series.
>
> Its purpose is to reduce the variance of the periodogram.
>
> > To smooth the spectrum, you need to either do ensemble averaging
> > (split the record into blocks, calculate the PSD of each and average),
> > or frequency averaging (accumulate the energy into large frequency
> > blocks), or both.
>
> If you do the maths, you might find that windowing in time
> domain is equivalent to convolution in frequency domain.
> Convolution is, in turn, equivalent to weighted averaging
> of the sequence, in this case the spectrum.
>
> Rune

Thanks to both of your replies. You are right H1 is Hl, sorry for the confusion. Secondly, I am more inclined to agree with Rune, window function does make the PSD curve smoother. But the question remains, how do I incorporate a cosine window function into my PSD program? I can do PSD plots with windowing manually in excel, but it taking me far too long. I was hoping Matlab would give me good results, but when I attempted it with the above code, it is giving me very scattered dominant frequencies. When I plotted the PSD curves, it was clear they weren't smooth. I know there are more exotic functions like Hamming, Blackman, Hann etc, but I want to start with a simple cosine window function above. Appreciate any assistance on this matter.

On Aug 29, 5:50 pm, "Safa " <enxs...@nottingham.ac.uk> wrote:
> I have a query concerning PSD. I am trying to use MATLAB for this, but the dominant frequencies I am getting are very scattered and I need to incorporate a windowing function to smooth out all the peaks. I would like to incorporate to begin with a simple cosine windowing function of the form w(kdt)=cos (Pi*kdt/2 *t). where dt is resolution of the time delay. I am unsure how to put this into my program, and this is simpler than the window function inbuilt in Matlab. The program is as follows (H1 in this >example has 30000 data points):
>
> Any suggestions to how to improve the program would be much appreciated.
> Many thanks, Safa


>
> H1=a(:,2);
> t=0.001:0.001:30;

t = (0:0.001:30-0.001)';

% See below:

N=length(t)
dt = t(2)-t(1)
T = N*dt % period of ifft(fft(H1))

t = (0:dt:T-dt)';
t = linspace(0,T-dt,N)';
t = dt*(0:N-1)';

> time_interval=t(2)-t(1);

Misleading notation. Notation implies t(end)-t(1)

> sizeHl=30000;

Misleading notation:

help size
help length

MYSTERY: When I cut and paste into my m-file
HI is mysteriously converted to H1...Why???

> Y = fft(Hl);
> fr = (1/time_interval)*(1:15000)/sizeHl;
> Pyy = Y.* conj(Y)/sizeHl;
> Pyy_interval=Pyy(2:15001);

Use of "interval' is misleading notation .

If you are not interested in the DC component,
first remove it in the time domain

H1 = H1-mean(H1);

Then you don't have to screw up your single sided
indexing by starting with "2"

> newcounter=1;
>
> [C,I] = max(Pyy_interval);

This will only find the first global maximum
To find all global maxima

tol = 1e-6 % set a tolerance
Imax = find(Pyy_interval > C-tol)

Need more complicated code to find all local peaks.

> f= (I/time_interval)/sizeHl;
> frequency =f'

Incorrect.

Fs = 1/dt
df = Fs/N
df = 1/T

frequency = df*(0:N/2)';
frequency = (0:df:Fs/2)';
frequency = linspace(0,Fs/2,N/2+1)';

freqmax = frequency(Imax)

To include windowing, replace H1 with w.*H1.

Hope this helps.

Greg

On Aug 29, 5:50 pm, "Safa " <enxs...@nottingham.ac.uk> wrote:
> I have a query concerning PSD. I am trying to use MATLAB for this, but the dominant frequencies I am getting are very scattered and I need to incorporate a windowing function to smooth out all the peaks. I would like to incorporate to begin with a simple cosine windowing function of the form w(kdt)=cos (Pi*kdt/2 *t). where dt is resolution of the time delay. I am unsure how to put this into my program, and this is simpler than the window function inbuilt in Matlab. The program is as follows (H1 in this >example has 30000 data points):
>
> Any suggestions to how to improve the program would be much appreciated.
> Many thanks, Safa
>
> H1=a(:,2);
> t=0.001:0.001:30;

t = (0:0.001:30-0.001)';

% See below:

N=length(t)
dt = t(2)-t(1) % time resolution
T = N*dt % period of ifft(fft(H1))

t = (0:dt:T-dt)';
t = linspace(0,T-dt,N)';
t = dt*(0:N-1)';

> time_interval=t(2)-t(1);

Misleading notation. Notation implies t(end)-t(1)

> sizeHl=30000;

Misleading notation:

help size
help length

MYSTERY: When I cut and paste into my m-file
HI is mysteriously converted to H1...Why???

> Y = fft(Hl);
> fr = (1/time_interval)*(1:15000)/sizeHl;
> Pyy = Y.* conj(Y)/sizeHl;
> Pyy_interval=Pyy(2:15001);

Use of "interval' is misleading notation .

If you are not interested in the DC component,
first remove it in the time domain

H1 = H1-mean(H1);

Then you don't have to screw up your single sided
indexing by starting with "2"

> newcounter=1;
>
> [C,I] = max(Pyy_interval);

This will only find the first global maximum
To find all global maxima

tol = 1e-6 % set a tolerance
Imax = find(Pyy_interval > C-tol)

Need more complicated code to find all local peaks.

> f= (I/time_interval)/sizeHl;
> frequency =f'

Incorrect.

Fs = 1/dt % Sampling frequency
df = Fs/N % Frequency resolution
df = 1/T

frequency = df*(0:N/2)';
frequency = (0:df:Fs/2)';
frequency = linspace(0,Fs/2,N/2+1)';

freqmax = frequency(Imax)

To include windowing, replace H1 with w.*H1.

Hope this helps.

Greg

Thanks very much for your comments, although there are what appears two identical replies? If I take on board your suggestions, is the new program as per below? It was a little difficult to follow, as there are tabs everywhere, so I tried to tidy it up. One query about the windowing equation. Appreciate your help.

H1=a(:,2); % Without windowing
%H1= w.*H1 % With windowing???
t = (0:0.001:30-0.001)';
N=length(t);
dt = t(2)-t(1); % time resolution
T = N*dt; % period of ifft(fft(H1))
t = (0:dt:T-dt)';
t = linspace(0,T-dt,N)';
t = dt*(0:N-1)';
S=size(H1);
Y = fft(H1);
fr = (1/dt)*(1:15000)/S;
Pyy = Y.* conj(Y)/S;
Pyy_interval=Pyy(2:15001);
H1 = H1-mean(H1);
newcounter=1;
tol= 1e-6 % set a tolerance
Imax = find(Pyy_interval > C-tol)
Fs = 1/dt % Sampling frequency
df = Fs/N % Frequency resolution
df = 1/T
frequency = df*(0:N/2)';
frequency = (0:df:Fs/2)';
frequency = linspace(0,Fs/2,N/2+1)';
freqmax = frequency(Imax);

On 10-08-30 11:22 AM, Greg Heath wrote:
> On Aug 29, 5:50 pm, "Safa "<enxs...@nottingham.ac.uk> wrote:

>> H1=a(:,2);

>> sizeHl=30000;

> MYSTERY: When I cut and paste into my m-file
> HI is mysteriously converted to H1...Why???

When I check the message here (via Usenet), the two-character variable name I
find ends in the digit one, not the letter capital-I .

I do not find any H followed by capital-I in what Safa posted, but I do see a
number of occurrences of H followed by lower-case L...

>> Y = fft(Hl);

including there. That variable, Hl, is not initialized.

On Aug 30, 1:31 pm, "Safa " <enxs...@nottingham.ac.uk> wrote:
> Thanks very much for your comments, although there are what
appears two identical replies? If I take on board your
suggestions, is the new program as per below? It was a little
difficult to follow, as there are tabs everywhere, so I tried
to tidy it up. One query about the windowing equation.
>Appreciate your help.
>
> H1=a(:,2); % Without windowing
> %H1= w.*H1 % With windowing???

Yes

> t = (0:0.001:30-0.001)';
> N=length(t);

N = length(t)

When coding and debugging I leave the
trailing semicolons off of scalar and small
matrix assignments so that they will print
to the command window.

> dt = t(2)-t(1); % time resolution
> T = N*dt; % period of ifft(fft(H1))
> t = (0:dt:T-dt)';
> t = linspace(0,T-dt,N)';
> t = dt*(0:N-1)';
> S=size(H1);
> Y = fft(H1);
> fr = (1/dt)*(1:15000)/S;

No.

To understand why, compare the results of
(No semicolons)

S = size(H1)
N = length(H1)

Also, the first frequency is zero.

Just modify the frequency statements below.

> Pyy = Y.* conj(Y)/S;
> Pyy_interval=Pyy(2:15001);
> H1 = H1-mean(H1);

No.

Remove mean before transforming

> newcounter=1;
> tol= 1e-6 % set a tolerance
> Imax = find(Pyy_interval > C-tol)

Where are Imax and C defined?

> Fs = 1/dt % Sampling frequency
> df = Fs/N % Frequency resolution
> df = 1/T
> frequency = df*(0:N/2)';
> frequency = (0:df:Fs/2)';
> frequency = linspace(0,Fs/2,N/2+1)';
> freqmax = frequency(Imax);

Please do not waste my time by posting code
that has not been tested!

If you have trouble getting it to run, include
the relevant error messages.

Again, if you remove the ending semicolons of
selected commands you can monitor the output.

Greg Heath

Sorry I was away on a conference and had to leave this for a few weeks. I am trying to follow your latest comments, and I am even more confused than before. You seem also to be contradicting your own code e.g. H1 = H1-mean(H1); One minute you say put it in, and next minute take it out. Please can I ask you for a HUGE favour. Could you kindly write your whole suggested program in one paragraph with no comments, and I will try it on my data? I will report back here on the results I get. I have an Excel macro that I can compare it with. Thanks for your assistance.

Here's the code as I understand it...
H1=a(:,2); % Without windowing
%w=cos (Pi*dt/2 *t) cosine window function
%H1= w.*H1 % With windowing
t = (0:0.001:30-0.001)';
N=length(t);
dt = t(2)-t(1); % time resolution
T = N*dt; % period of ifft(fft(H1))
t = (0:dt:T-dt)';
t = linspace(0,T-dt,N)';
t = dt*(0:N-1)';
N1=length(H1);
Y = fft(H1);
fr = (1/dt)*(1:15000)/N1;
Pyy = Y.* conj(Y)/N1;
Pyy_interval=Pyy(2:15001);
H1 = H1-mean(H1);
newcounter=1;
[C,I] = max(Pyy_interval);
tol= 1e-6 % set a tolerance
Imax = find(Pyy_interval > C-tol)
Fs = 1/dt % Sampling frequency
df = Fs/N % Frequency resolution
df = 1/T
frequency = df*(0:N/2)';
frequency = (0:df:Fs/2)';
frequency = linspace(0,Fs/2,N/2+1)';
freqmax = frequency(Imax);

"Safa " <enxss10@nottingham.ac.uk> wrote in message <i80f14$3ts$1@fred.mathworks.com>...
> Sorry I was away on a conference and had to leave this for a few weeks. I am trying to follow your latest comments, and I am even more confused than before. You seem also to be contradicting your own code e.g. H1 = H1-mean(H1); One minute you say put it in, and next minute take it out. Please can I ask you for a HUGE favour. Could you kindly write your whole suggested program in one paragraph with no comments, and I will try it on my data? I will report back here on the results I get. I have an Excel macro that I can compare it with. Thanks for your assistance.
>
> Here's the code as I understand it...
> H1=a(:,2); % Without windowing
> %w=cos (Pi*dt/2 *t) cosine window function
> %H1= w.*H1 % With windowing
> t = (0:0.001:30-0.001)';
> N=length(t);
> dt = t(2)-t(1); % time resolution
> T = N*dt; % period of ifft(fft(H1))
> t = (0:dt:T-dt)';
> t = linspace(0,T-dt,N)';
> t = dt*(0:N-1)';
> N1=length(H1);
> Y = fft(H1);
> fr = (1/dt)*(1:15000)/N1;
> Pyy = Y.* conj(Y)/N1;
> Pyy_interval=Pyy(2:15001);
> H1 = H1-mean(H1);
> newcounter=1;
> [C,I] = max(Pyy_interval);
> tol= 1e-6 % set a tolerance
> Imax = find(Pyy_interval > C-tol)
> Fs = 1/dt % Sampling frequency
> df = Fs/N % Frequency resolution
> df = 1/T
> frequency = df*(0:N/2)';
> frequency = (0:df:Fs/2)';
> frequency = linspace(0,Fs/2,N/2+1)';
> freqmax = frequency(Imax);

Hi Safa, Do you have the Signal Processing Toolbox? If so, you can use spectral analysis objects to get a PSD estimate very easily. I'm not sure if when you say "smooth" out the peaks you are referring to the variance of the periodogram estimate or something else. You can greatly reduced the variance by an averaging technique as Tideman mentioned in an earlier post in this thread.

Fs = 1000;
t = linspace(0,1,1000);
x = cos(2*pi*100*t)+0.5*sin(2*pi*200*t)+randn(size(t));
% Get periodogram estimate
hper = spectrum.periodogram;
psdper = psd(hper,x,'NFFT',length(x),'Fs',Fs);
plot(psdper)
% Get Welch estimate using Hamming windwo
 hw = spectrum.welch('Hamming',100);
psdwelch = psd(hw,x,'NFFT',length(x),'Fs',Fs);
figure;
plot(psdwelch)

Notice the difference in the variability of the PSD estimate compared with the Welch estimate.

Wayne

On Sep 29, 6:36 pm, "Safa " <enxs...@nottingham.ac.uk> wrote:
> Sorry I was away on a conference and had to leave this for a few weeks. I am trying to follow your latest comments, and I am even more confused than before. You seem also to be contradicting your own code e.g. H1 = H1-mean(H1); One minute you say put >it in, and next minute take it out.

Your sentences are not wrapping when pasted into Notepad.
Please make sure carriage returns are included.

It depends on what YOU want. Do you understand the difference?
If you don't, that is the question you should be asking
...or, just run a simple test case

T = 2*pi, N = 8, dt = T/N
t = dt*(0: N-1)'; H1 = 1 + sin(t);
Y1 = fft(H1)/N ; Y2 = fft(H1-mean(H1))/N;

Now what is the SIGNIFICANT difference between
Y1 and Y2 (and what relationship is that to
fft(ones(N,1)/N ?)

dY = Y1-Y2

>Please can I ask you for a HUGE favour. Could you kindly
> write your whole suggested program in one paragraph
> with no comments, and I will try it on my data?

My goal is to give you as many hints as you need
so that you can understand enough to be able to
write your own program

Suggestion for coding and debugging:
1. Use a small test value for N, e.g., 8 if original N
is even or 9 if original N is odd
2. Truncate corresponding variables; e.g., t=t(1:N)
and H1 = H1(1:N)
3. Remove the ending semicolons of selected
commands during coding and debugging, so that
you can easily monitor the calculations to make
sure they are not unreasonable.

> I will report back here on the results I get. I have an
> Excel macro that I can compare it with.
>Thanks for your assistance.
>
> Here's the code as I understand it...
> H1=a(:,2); % Without windowing
> %w=cos (Pi*dt/2 *t) cosine window function
> %H1= w.*H1 % With windowing
> t = (0:0.001:30-0.001)';
> N=length(t);

Is it equal to N1 = length(H1)?
Is it even or odd?

If it is the same, eliminate confusion by
using either N or N1, not both, in all of the
following

> dt = t(2)-t(1); % time resolution
> T = N*dt; % period of ifft(fft(H1))
> t = (0:dt:T-dt)';
> t = linspace(0,T-dt,N)';
> t = dt*(0:N-1)';

The last 3 statements are equivalent.
Choose which format you like.

> N1=length(H1);
> Y = fft(H1);

In general:

If you are going to ignore Y(1), it is better to remove
the mean from H1 so that the sinc function sidelobes
from Y(1)*ones(N,1) do not contaminate Y(2:N).

> fr = (1/dt)*(1:15000)/N1;

corresponds to the 3 equations for frequency
below with frequency(1) removed

> Pyy = Y.* conj(Y)/N1;
> Pyy_interval=Pyy(2:15001);

In general, the single-sided form of Pyy is

if mod(N,2)~=0 % N is odd
     Pyyss = [ Pyy(1) ; 2*Pyy(2:(N+1)/2) ];
else % N is even
     Pyyss = [ Pyy(1) ; 2*Pyy(2:N/2); Pyy(N/2+1)];
end

> H1 = H1-mean(H1);

see above comments

> newcounter=1;

What is being counted?

> [C,I] = max(Pyy_interval);
> tol= 1e-6 % set a tolerance
> Imax = find(Pyy_interval > C-tol)
> Fs = 1/dt % Sampling frequency
> df = Fs/N % Frequency resolution
> df = 1/T
> frequency = df*(0:N/2)';
> frequency = (0:df:Fs/2)';
> frequency = linspace(0,Fs/2,N/2+1)';

The last 3 statements are equivalent.
Choose which format you like.

> freqmax = frequency(Imax);


Hope this helps.

Greg

Wayne may I say that is a fantastic example you’ve given, and it illustrates precisely my problem and what I am trying to do!! Yes I have the signal processing toolbox, although I must admit I haven't used it. In your example, you can see how much smoother the PSD is when you apply a window function. The question how can I apply a simple cosine window function using your approach? Is there a way to define it manually? I have been trying to do this with Greg, and it has been a painful exercise. While I highly appreciate Greg’s wonderful efforts, it is getting into fairly advanced signal processing techniques. But I like Greg’s methodology to get me thinking and I will keep trying, until I get my head round everything he has suggested. I have a lot of data to process and it is a nightmare to do this in Excel. I need a more automated process this is why I turned to
Matlab. I really appreciate everyone's help. Many thanks for your assistance.

"Safa " <enxss10@nottingham.ac.uk> wrote in message <i83aqr$e7r$1@fred.mathworks.com>...
> Wayne may I say that is a fantastic example you’ve given, and it illustrates precisely my problem and what I am trying to do!! Yes I have the signal processing toolbox, although I must admit I haven't used it. In your example, you can see how much smoother the PSD is when you apply a window function. The question how can I apply a simple cosine window function using your approach? Is there a way to define it manually? I have been trying to do this with Greg, and it has been a painful exercise. While I highly appreciate Greg’s wonderful efforts, it is getting into fairly advanced signal processing techniques. But I like Greg’s methodology to get me thinking and I will keep trying, until I get my head round everything he has suggested. I have a lot of data to process and it is a nightmare to do this in Excel. I need a more automated process this is why I turned to
> Matlab. I really appreciate everyone's help. Many thanks for your assistance.

Hi Safa, Yes, Greg certainly invests a lot of his time and experience to help people with questions.

If you wish to apply just a single window by using a spectral analysis object, you can do that with spectrum.periodogram. However, the reduction in the variance of the PSD estimate will not match that of the Welch estimator. To use a single window:

Fs = 1000;
t = linspace(0,1,1000);
x = cos(2*pi*100*t)+0.5*sin(2*pi*200*t)+randn(size(t));
% Using Hamming window
h = spectrum.periodogram('Hamming');
psdhamm = psd(h,x,'NFFT',length(x),'Fs',Fs);
plot(psdhamm)

There are a number of windows supported in the Signal Processing Toolbox. To see all the possible built-in windows, you can enter:

set(spectrum.periodogram,'WindowName')

Note the last option. If you really wish to use your own window function, that is also possible by using the 'user defined' option. If that is the case, write back and I'll show you how to do that. I suspect that for your purposes, one of the windows on that extensive list will suffice.

Wayne

Wayne, Thank you very much for your help and I agree with you, the Welch spectrum is really powerful, I am amazed how it is filtering the PSD. If I want to automate your program and run say 10 data sets, how can I get Matlab to plot all 10 PSDs on one sheet (separate graphs), and to also generate a report giving me the maximum frequency for each PSD? For the code I wrote below, in the workspace, it only saves data for the last run. I've worked out how to do the PSD manually using the signal processing toolbox, which is a really useful facility to have and having played with all the windows, I found Hanning would be a good choice. I’m working with a file that contains both numbers and words, and matlab does not like that at all so need to filter the data. If you have time, I would be interested in the “user defined” command, for a cosine window of the form w(kdt)=cos
(Pi*kdt/2 *t). Appreciate your help.


for run=002:012
    
if run<=9
fname=['G:\run00',num2str(run),'.epst'];
elseif 10<=run<=99
fname=['G:\run0',num2str(run),'.epst'];
end

[h, a] = mhdrload(fname);
H=a(:,2);
H1=H/100;
 
Fs = 1000;
% Get Welch estimate using Hamming window
 hw = spectrum.welch('Hann',100);
psdwelch = psd(hw,H1,'NFFT',length(H1),'Fs',Fs);
plot(psdwelch)

end

"Safa " <enxss10@nottingham.ac.uk> wrote in message <i8ak18$n09$1@fred.mathworks.com>...
> Wayne, Thank you very much for your help and I agree with you, the Welch spectrum is really powerful, I am amazed how it is filtering the PSD. If I want to automate your program and run say 10 data sets, how can I get Matlab to plot all 10 PSDs on one sheet (separate graphs), and to also generate a report giving me the maximum frequency for each PSD? For the code I wrote below, in the workspace, it only saves data for the last run. I've worked out how to do the PSD manually using the signal processing toolbox, which is a really useful facility to have and having played with all the windows, I found Hanning would be a good choice. I’m working with a file that contains both numbers and words, and matlab does not like that at all so need to filter the data. If you have time, I would be interested in the “user defined” command, for a cosine window of the form w(kdt)=cos
> (Pi*kdt/2 *t). Appreciate your help.
>
>
> for run=002:012
>
> if run<=9
> fname=['G:\run00',num2str(run),'.epst'];
> elseif 10<=run<=99
> fname=['G:\run0',num2str(run),'.epst'];
> end
>
> [h, a] = mhdrload(fname);
> H=a(:,2);
> H1=H/100;
>
> Fs = 1000;
> % Get Welch estimate using Hamming window
> hw = spectrum.welch('Hann',100);
> psdwelch = psd(hw,H1,'NFFT',length(H1),'Fs',Fs);
> plot(psdwelch)
>
> end


Hi Safa, if you look at the contents of psdwelch, just enter

>>psdwelch

You will see that psdwelch contains fields called Data and Frequencies among others.

You can easily extract those fields and save them. For example, the following code saves the PSD estimates for 10 time series in a matrix, identifies the maximum PSD value for each PSD estimate, and finds the frequency corresponding to the maximum. The frequency corresponding to the maximum is stored in the variable maxfreq.

Fs = 1000;
t = linspace(0,1,1000);
psddata = zeros(501,10);
h = spectrum.welch;
for J=1:10
x = cos(2*pi*100*t)+randn(size(t));
psdwelch = psd(h,x,'NFFT',length(x),'Fs',Fs);
psddata(:,J) = psdwelch.Data;
[Y,I] = max(psdwelch.Data);
maxfreq(J) = psdwelch.Frequencies(I);
end

Wayne

"Safa " <enxss10@nottingham.ac.uk> wrote in message <i5ekmr$4u$1@fred.mathworks.com>...
> I have a query concerning PSD. I am trying to use MATLAB for this, but the dominant frequencies I am getting are very scattered and I need to incorporate a windowing function to smooth out all the peaks. I would like to incorporate to begin with a simple cosine windowing function of the form w(kdt)=cos (Pi*kdt/2 *t). where dt is resolution of the time delay. I am unsure how to put this into my program, and this is simpler than the window function inbuilt in Matlab. The program is as follows (H1 in this example has 30000 data points):
>
> H1=a(:,2);
> t=0.001:0.001:30;
> time_interval=t(2)-t(1);
> sizeHl=30000;
> Y = fft(Hl);
> fr = (1/time_interval)*(1:15000)/sizeHl;
> Pyy = Y.* conj(Y)/sizeHl;
> Pyy_interval=Pyy(2:15001);
> newcounter=1;
>
> [C,I] = max(Pyy_interval);
> f= (I/time_interval)/sizeHl;
>
> frequency =f'
>
> Any suggestions to how to improve the above program would be much appreciated.
> Many thanks, Safa

help tukeywin

Ok I am baffled again, I tried using the welch method, and it keeps giving me the same maximum frequency for all different runs. It is not picking up local maximum frequencies. So I have gone back to my initial program, and incorporated what Greg has suggested. It works. However, I have one final query about the maximum frequency it reads. For one run with 30000 data points (30 seconds of data @ Fs=1000Hz), the matrix for frequency = (0:df:Fs/2)'; is 15001xdouble and starts at zero. The matrix for fr = (1/dt)*(1:15000)/N is 15000xdouble and doesn't start at zero. The matrix for Pyy interval is 15000xdouble. In the trial run, it says it is picking out row 29 as maximum frequency. Is it reading the same value in the same row for both Pyy interval and frequency matrices? Why are the two matrices not the same size and do we need to adjust them to make them same size? Appreciate your help
once more.

Ok I will make my question simpler. I've tidied up the code as follows:

Hl=a(:,2);
t = (0:0.001:30-0.001)';
N=length(t);
dt = t(2)-t(1);
T = N*dt; % period of ifft(fft(H1))
t = (0:dt:T-dt)';
%t = linspace(0,T-dt,N)';
%t = dt*(0:N-1)';
N=length(Hl);
Y = fft(Hl);
fr = (1/dt)*(1:15000)/N;
Pyy = Y.* conj(Y)/N;
Pyy_interval=Pyy(2:15001);
[C,I] = max(Pyy_interval);
tol= 1e-6; % set a tolerance
Imax = find(Pyy_interval > C-tol);
Fs = 1/dt; % Sampling frequency
df = Fs/N; % Frequency resolution
df = 1/T;
frequency = df*(0:N/2)';
%frequency = (0:df:Fs/2)';
%frequency = linspace(0,Fs/2,N/2+1)';

My question concerns the final part of the code which is the maximum frequency. Greg suggested
freqmax = frequency(Imax)

Should it not be
freqmax=fr (Imax)

If not why not? The size of matrix of "fr" corresponds with the size of matrix of "Pyy_interval" which 1x15000. However "frequency" generates 1x15001 matrix and starts at zero. I hope my query is now clearer and thank you again so much for your help.

I apologise slightly corrected code (I must ignore Y(1), it is so much higher than the other values, so taken on board your suggestion to take out mean value for Hl)

Hl=a(:,2); % Without windowing
%w=cos (Pi*dt/2 *t) cosine window function
%Hl=w.*H1 % With windowing
Hl=Hl-mean(Hl);
t = (0:0.001:30-0.001)';
N=length(t);
dt = t(2)-t(1);
T = N*dt; % period of ifft(fft(H1))
t = (0:dt:T-dt)';
%t = linspace(0,T-dt,N)';
%t = dt*(0:N-1)';
N=length(Hl);
Y = fft(Hl);
fr = (1/dt)*(1:15000)/N;
Pyy = Y.* conj(Y)/N;
Pyy_symmetry=Pyy(2:15001);
[C,I] = max(Pyy_symmetry);
tol= 1e-6; % set a tolerance
Imax = find(Pyy_symmetry > C-tol);
Fs = 1/dt; % Sampling frequency
df = Fs/N; % Frequency resolution
df = 1/T;
%frequency = df*(0:N/2)';
frequency = (0:df:Fs/2)';
%frequency = linspace(0,Fs/2,N/2+1)';
freqmax = frequency(Imax) OR freqmax=fr(Imax)????

Would really appreciate some help with the above program. I think the bulk of the program is now OK, I just need clarification to a couple of lines in the code. Thanks in anticipation!

On Oct 5, 8:20 pm, "Safa " <enxs...@nottingham.ac.uk> wrote:
> I apologise slightly corrected code (I must ignore Y(1), it is so much higher than the other values, so taken on board your suggestion to take out mean value for Hl)
>
> Hl=a(:,2); % Without windowing
> %w=cos (Pi*dt/2 *t) cosine window function
> %Hl=w.*H1 % With windowing
> Hl=Hl-mean(Hl);
> t = (0:0.001:30-0.001)';

It is better to use symbols than numbers.
See the expression for t 4 equations below.

> N=length(t);

How do you know that t and HI have the same length?

> dt = t(2)-t(1);
> T = N*dt; % period of ifft(fft(H1))
> t = (0:dt:T-dt)';

You already defined t above

> %t = linspace(0,T-dt,N)';
> %t = dt*(0:N-1)';
> N=length(Hl);
> Y = fft(Hl);
> fr = (1/dt)*(1:15000)/N;

What combination of parameters is equal to 15000?

> Pyy = Y.* conj(Y)/N;

Pyy = Y.*conj(Y)/(N*Fs);

> Pyy_symmetry=Pyy(2:15001);

Lousy notation. Use Pyyss for single-sided
Double all components except the last to account
fot the negative freq

> [C,I] = max(Pyy_symmetry);
> tol= 1e-6; % set a tolerance
> Imax = find(Pyy_symmetry > C-tol);
> Fs = 1/dt; % Sampling frequency
> df = Fs/N; % Frequency resolution
> df = 1/T;
> %frequency = df*(0:N/2)';
> frequency = (0:df:Fs/2)';
> %frequency = linspace(0,Fs/2,N/2+1)';
> freqmax = frequency(Imax) OR freqmax=fr(Imax)????

On Oct 18, 12:10 am, Greg Heath <he...@alumni.brown.edu> wrote:
> On Oct 5, 8:20 pm, "Safa " <enxs...@nottingham.ac.uk> wrote:
>
> > I apologise slightly corrected code (I must ignore Y(1), it is so much higher than the other values, so taken on board your suggestion to take out mean value for Hl)
>
> > Hl=a(:,2); % Without windowing
> > %w=cos (Pi*dt/2 *t) cosine window function
> > %Hl=w.*H1 % With windowing
> > Hl=Hl-mean(Hl);
> > t = (0:0.001:30-0.001)';
>
> It is better to use symbols than numbers.
> See the expression for t 4 equations below.
>
> > N=length(t);
>
> How do you know that t and HI have the same length?
>
> > dt = t(2)-t(1);
> > T = N*dt; % period of ifft(fft(H1))
> > t = (0:dt:T-dt)';
>
> You already defined t above
>
> > %t = linspace(0,T-dt,N)';
> > %t = dt*(0:N-1)';
> > N=length(Hl);
> > Y = fft(Hl);
> > fr = (1/dt)*(1:15000)/N;
>
> What combination of parameters is equal to 15000?
>
> > Pyy = Y.* conj(Y)/N;
>
> Pyy = Y.*conj(Y)/(N*Fs);
>
> > Pyy_symmetry=Pyy(2:15001);
>
> Lousy notation. Use Pyyss for single-sided
> Double all components except the last to account
> fot the negative frequency power..
>
> > [C,I] = max(Pyy_symmetry);
> > tol= 1e-6; % set a tolerance
> > Imax = find(Pyy_symmetry > C-tol);
> > Fs = 1/dt; % Sampling frequency
> > df = Fs/N; % Frequency resolution
> > df = 1/T;
> > %frequency = df*(0:N/2)';
> > frequency = (0:df:Fs/2)';
> > %frequency = linspace(0,Fs/2,N/2+1)';

The last 6 statements were given to
help explain the definition of the
single-sided frequency grid. Modified
to exclude zero frequency, the 5th one
would replace your uninformative equation
for fr (What does the second letter "r" denote?)

> freqmax = frequency(Imax) OR freqmax=fr(Imax)????

Since you excluded the DC component in the
definition of Pyy_symmetry, it is the latter.

Hope this helps.

Greg