How to preprocess signals prior cross-correlation? Problem Description

I have two devices which are measuring the same inputs and are using them to calculate the same (non-periodic) output signal. The process generating the inputs is deterministic.
Both devices are sampling with different rates (each > 5 kHz).
The capture duration of one device may be several seconds longer than the other one. Total capture duration is > 60s.
Since the calculation of the output signals is done with different algorithms they are looking quite different at first sight but they do look very similar after applying some sort of linear low-pass filter.

My goal is to find out the time delay between those two output signals x and y.
Proposal for Solution

Crop the longer signal such that x and y have the same capture duration.
Resample x and y such that they have the same number of samples.
Filter x and y with the same linear low-pass filter.
Remove linear trends by applying MATLAB's detrend function to x and y.
Apply MATLAB's finddelay on x and y.

Questions

May cropping the longer signal corrupt the result of the cross-correlation?

What is the best approach for resampling x and y?

As I mentioned, x and y are only looking very similar
after filtering. Is filtering signals prior cross-correlation
advisable? If yes, are there recommended filter types for that
purpose?

Is it a must to detrend data prior cross-correlation or might it be enough to remove the mean? In my case there may be multiple linear trends and it could be hard to find all breakpoints as depicted here.

Desired Time Delay between two (filtered) Output Signals

What you are describing is called (/ has been solved via) "Dynamic Time Warping".

In DTW, you calculate a series of simple operations between two sequences that if applied on one sequence they would "bring it" directly "over" the other. Of course, there might be more than one ways to do this in which case you are selecting the one with the least cost. This "cost of distortion" is now your similarity.

For more information on DTW, please see this, this and this. More relevant to your application might be this publication. For a practical example, please see this link.

With regards to some of the rest of your questions:

  May cropping the longer signal corrupt the result of the cross-correlation?

Yes, it would...

  What is the best approach for resampling x and y?

...because both cropping and re-sampling assume a linear speed up. This might not be the case across the entire length of your signal.

  As I mentioned, x and y are only looking very similar after filtering. Is filtering signals prior cross-correlation advisable?

It depends a lot on the application. However, filtering two signals $x,y$ with the same impulse response $h$ is expected to increase their cross correlation because after convolution, both $x,y$ time instances would now have the same $h$ imprint on them.

EDIT: Following a brief review of the data provided with the question, it seems that these are some sort of counts (?) or other digital / impulsive data. Therefore, the question "...is filtering signals prior to cross-correlation advisable?" is ill posed. The fact remains that if you filter two signals through the same filter, their cross-correlation will increase. But here, it seems that the signals do not exist in their final form yet. To clarify: We have to agree on which one is the signal. If you were dealing (for example) with a pulse width modulation (PWM) derived sound signal, we would not try to correlate the two PWM pulse trains that produce it but their integrated versions because THAT is what we consider as our signal. Consequently, in this particular case we are dealing with here, we have to understand what is considered as "signal". The only thing we know about it is that it comes out of some process. That's fine but, at the same time, if you have to integrate the data to PRODUCE the signal of interest then the process itself or its context would tell you what sort of integration you have to apply. It might not even be called "low pass filtering". In conclusion: If your signal has to be integrated first, then don't try to compare the "pulse sequences". To get back to the PWM-Audio case, it is like asking "Should I integrate the PWM pulses before cross correlation or would this spoil my signal?". Of course you have to apply the integration because otherwise, you don't have the signal to work with. Decide on the parameters of integration, apply them to your original signal and then, if you still have drifts and compressions and dilations of features in the signal, see if you can use DWT to establish the time delay.

  Is it a must to detrend data prior cross-correlation or might it be enough to remove the mean?

Again, depends on the application and the amount of drift expected in the data but in general, it is a good idea to remove any long term drifts to avoid the cross-correlation being thrown off by the fact that the trends seem to match despite those little "details" on top of them. Where "details", of course, might be your desired signal. This is illustrative but a hypothetical scenario, your data might not even include such a case, it simply shows how detrending might help.

Hope this helps.

EDIT (2):

In the light of more information, I have to say that this does not look like a task for DTW neither cross-correlation the way it is described.

What is attempted here is an estimation of the delay between two waveforms by aligning "landmark" points. NOT the waveforms themselves. Personally, I am struggling to see why the marked points in the signal should align at all. Therefore, what I would suggest is this:

Pick the higher frequency waveform
Mark the peak you are after
Extract $\pm m$ points (or ms) around it. For example, if the peak is at sample 4500, extract x[4500-200:4500+200]. Call this your template $q$.
Subtract the mean of $q$ from $q$
Downsample $q$ so that it matches the sampling rate of the lower frequency signal. Call this $q_d$.

This is now your template. It is shorter than the signal you are looking for and you can "slide" it over the lower rate signal and estimate the cross correlation. At the point that the cross correlation peaks is where this second waveform is, in the lower rate signal. Prior to applying cross correlation, you might want to be subtracting the mean of the signal portion you are about to compare with your template or instead of removing the mean, you might want to be applying detrending. Because the duration of the signal is small, simple linear detrending might work good enough.

This is an attempt to align specific points of interest in the signal. Not the waveform as a whole

Hope this helps.

EDIT (3):
Ideally, a large amount of the information provided in the comments here, should appear in the original question and then perhaps we can delete all of these comments.

If you can control the load of the inverter, maybe you can insert a "sync event". Apply a known load on and off after the inverter has reached its nominal operation to create a "click" with a predictable response. You can use that "click" later to sync the two signals. If you have this capability, then you can average the responses from a number of "clicks" so that you have a good idea of what you are looking for, via cross correlation, as the "sync pulse".

Of course, this covers what happens when the two signals are supposed to have a constant time invariant delay.


Interleaving or randomisation? I am trying to implement a digital communication system as a learning exercise. For error correction, I am using two concatenated codes (RS and convolutional). However, I read that, in order to reduce the damage caused by burst errors (and in some cases controlling the DC level of the output signal), it is a good practice to add an interleaver/scrambler between them.

I know how an interleaver and a scrambler work, but I do not get what is the point of choosing one or another (if they are not both required). In what scenarios would one work better than the other and why?

So, the technical difference between a scrambler and an interleaver:

The interleaver is just a permutation of input symbols; a resorting.
The scrambler is an operation designed to modify the transmitted symbol in such a way that they are not present in the output.

Thus, "scrambler" is a way less precise defined term. I usually define it as "operation that is easy to invert if you know the transmit operation, and typically fulfills more of a statistical purpose than a cryptographic one".

The typical purpose of a scrambler is to break up periodicities in the original symbol stream; this can, for example, be done by combining with a pseudorandom sequence. To de-scramble, the receiver would do the inverse combination (in case of XOR, it's autoinverse) with the same sequence. Often, that sequence is longer than a "block" of transmitted signals.

That way, you're basically applying an uncorrelated signal to your data signal; even if you had spikes in the spectrum of your data signal, because it had repetitiveness in it, these will (in expectation) be gone after scrambling, and you get a nice, white spectrum.

A white spectrum is important:

you will need to adhere to legal restrictions on power, and these usually specify "maximum power within a small bandwidth". So, if you want to use the most of your expensive spectrum, you want to have a nice PSD without peaks.
from frequency-selective channel point of view, having peaks in the spectrum is a dangerous game: if a deep fade happens to fall on the peak, you lose a lot of signal. You can "distribute" your risk better if you evenly use the channel.
from a synchronization point of view, for many synchronization algorithms, it's important that

the spectrum is symmetric (e.g. for the band-edge frequency locked loop (FLL)) or
that all symbols are equally likely (you'd have a peak at 0 Hz in your spectrum if you have an "average != 0" signal) (important for PLLs, otherwise the assumption that all symbols are equally likely breaks, and that breaks the optimality of "correcting" towards the closest correct phase)
that there's no runs of constant symbols, even if that happens likely in the input data (e.g. constant-rate digital broadcasting: when there's no bits to transmit, because e.g. everything is silent or the video codec only needs to compress a still image, then you have to pad with 0s, so that you keep your symbol rate) (extremely crucial for symbol timing recovery to have enough symbol transitions)

The interleaver, on the other hand, doesn't change the fact that your data signal might have periodic components: sure, you're moving them to different positions, but they're still there; you just "moved" the spikes in the spectrum in phase, you have little chance of actually getting rid of them.

On the other hands, if you have a burst error on say symbols Nr. 100 to Nr. 200, your XOR-scrambler would still have that burst error in the same range. After a sufficiently long interleaver, the same burst error would be distributed to more or less isolated singular errors, which might be much better for the channel  decoder down the processing line.

So, all in all:

interleaver against burst errors
scrambler against unwanted regularities in the signal.


Finding inverse Z Transform Find the inverse Z transform:

I have done the solution but my answer does not match with the one given in the textbook. What I may have done wrong?

HINT: Your partial fraction expansion is wrong; you should always check it by recombining the factors. The correct expansion is

$$X(z)=\frac{z^{-1}-\frac12}{\left(1-\frac12 z^{-1}\right)^2}=\frac{\frac32}{\left(1-\frac12 z^{-1}\right)^2}-\frac{2}{\left(1-\frac12 z^{-1}\right)}\tag{1}$$

You should try to derive this result yourself.

BUT: note that you don't even need partial fractions here because if you know (or look up) the $\mathcal{Z}$-transform correspondence

$$\frac{az^{-1}}{\left(1-a z^{-1}\right)^2}\Longleftrightarrow na^nu[n]\tag{2}$$

then you can directly apply it by writing $X(z)$ as

$$X(z)=2\frac{\frac12z^{-1}}{{\left(1-\frac12 z^{-1}\right)^2}}-z\frac{\frac12z^{-1}}{{\left(1-\frac12 z^{-1}\right)^2}}\tag{3}$$

Finally, note that the given result is wrong. It should be

$$x[n]=2n\left(\frac12\right)^nu[n]-(n+1)\left(\frac12\right)^{n+1}u[n+1]\tag{4}$$


White noise gain of uniform spline interpolation Uniform spline interpolation is fully defined by its continuous-domain impulse response. I noticed that the integral of the square of the impulse response over its support gives the average amplification of the power of white noise error in the samples.

What is this figure or its square root called? Is "white noise gain" $(\operatorname{WNG})$ appropriate for the square root?

For example for linear interpolation:

Figure 1. Linear interpolation impulse response

$$\operatorname{WNG} = \sqrt{\int_{-1}^1\left(\left\{\begin{array}x+1&\text{if }x<0\\1-x&\text{otherwise}\end{array}\right.\right)^2dx} = \sqrt{\frac{2}{3}}\approx0.8164965809$$

This figure is the squared norm or energy of the impulse response:

$$
\|h(t)\|^2=\int_{\mathrm{R}} |h(t)|^2 dt
$$

As you point out correctly, this is the noise amplification. Considering a signal

$$
r(t) = x(t) + n(t)
$$
with $x(t)$ is the signal of interest and $n(t)$ is AWGN with variance $\sigma^2$. Then, after filtering

$$
\begin{align}
y(t) &= h(t)*r(t)\\
&=h(t)*x(t) + \int_{\mathrm{R}} h(t-\tau)n(\tau)d\tau\\
&=h(t)*x(t) + \tilde{n}(t)
\end{align}
$$

with $\tilde{n}(t)$ being the filtered noise. Now, considering

$$
\begin{align}
E[\tilde{n}(t)\tilde{n}^*(t)] &= \int_\mathrm{R}\int_\mathrm{R} h(t-\tau_1)h^*(t-\tau_2)E[n(\tau_1)n(\tau_2)]d\tau_1 d\tau_2\\
&=\sigma^2\int_{\mathrm{R}^2}h(t-\tau_1)h^*(t-\tau_2)\delta(\tau_1-\tau_2)d\tau_1d\tau_2\\
&=\sigma^2\int_\mathrm{R}h(t-\tau_1)h^*(t-\tau_1)d\tau_1\\
&=\sigma^2\|h(t)\|^2
\end{align}
$$

where we have used that the noise is white, i.e. its autocorrelation is a delta-function. Note that after filtering generally the noise is not white anymore, but colored.


Could multiple echo reflections be combined to improve signal to noise ratio? I am working on a system processing ultrasonic echo reflection. The transmitter creates a short ultrasonic pulse and the receiver analyses the reflection.

If I send two pulses (back to back or with some small defined delay) I receive two superimposed reflections.
I can process these superimposed reflection with the IIR filter described in this post:

Trivial echo separation

to restore the single reflection.

My question is, could I get a signal with better signal to noise ratio by using two or more transmit pulses compared to the signal from a single pulse?

You need to be careful to look at your bunch of impulses in the fequency domain, to make sure you are really transmitting more energy at all frequencies. Your chain of 3 identical impulses will contain zero energy at certain frequencies. As a result, your echo-separation method (Tk = Sk – (Tk-10) – (Tk-20)) will have infinite gain at those frequencies (this is an IIR filter with 2 poles on the unit-circle), and hence it will be unstable.

As an example, I built 4 different "pulse train" signals (in Matlab), by the following method (an I've plotted the magnitude of the FFT):

>> S1 = [1;0;0;0;0;0;0;0;0;0;0];
>> S2 = [1;0;0;0;0;1;0;0;0;0;0];
>> S3 = [1;0;0;0;0;1;0;0;0;0;1];
>> S4 = [1;0;0;0;0;1;0;0;0;0;-1];
>> plot( abs([fft(S1,256),fft(S2,256),fft(S3,256),fft(S4,256)]) );

Note that S1 has a constant magnitude respnnse (equal to 1), as expected. S2 (made up of 2 impulses) has a magnitude response that goes to 0 at several frequencies (even though, on average, it has twice the power of S1).

S3 has three times the average power as S1, but unfortunately it also has zero energy at some frequencies.

S4 is built with a -1 in the last sample. As a result, it still has 3 times the average energy compared to S1, with the added property that is it never less than unity energy in any one frequency band.

The (simplistic) recovery-filter, for use with S4, would be :  Tk = Sk – (Tk-5) + (Tk-10)), but unfortunately this is also an unstable filter (having some poles outside the unit circle). The trick to making a filter like this stable is split it into two parts.

First, you run a regular IIR filter over the signal as follows:

$$T'_k = S_k + g \times T'_{k-5}$$

where $g=\frac{\sqrt5-1}{2}$. Next, you run this filter:

$$T_k = T'_k - g \times T_{k+5}$$

and the trick, with this second filter, is to compute the samples of $T_k$ in reverse order (if you want to compute 100 samples, start with sample 99, and work your way down to sample 0).

The overarching principles here are:

Spread your energy out over multiple samples, in the transmitted signal, to allow you to send more energy
Try to always ensure there is lots of energy in all frequencies of the transmitted signal (look at the FFT of the signal)
Build an inverse filter

My inverse filter method shown here is a bit perverse, using a time-reversed IIR filter. If you have plenty of CPU cycles, you can compute the correct inverse filter by building it as a long FIR filter (in theory, you need an infinitely long filter, but a long enough FIR will approximate it OK).


Simplified equations for FM (frequency modulation) synthesis and PM (phase modulation)? Clicking on the links I find on the web it seems many authors don't express them
neatly enough for me or they're even confusing the two. Can someone please clear this up?

Is frequency modulation

\begin{align}
&\rm  \sin((frequency+modulator)time-phase)\\
&\qquad\quad\textrm{or is it}\\
& \rm  \sin(frequency*time+modulator-phase)\quad ?
\end{align}

Ugh!

If $t$ is time, $s(t)$ is the (appropriately scaled) signal, $\omega_0$ is the angular frequency, and $\phi_0$ is a phase offset, then phase modulation is
$$ t \mapsto \sin(\omega_0 t + \phi_0 + s(t)) $$
Different signs for $s(t)$ and/or $\phi_0$ may be used, depending on conventions and context.

Frequency modulation is more complex to write down, because we want $s(t)$ to vary the derivative of the argument to the sine as a function of time. We get something like:
$$ t \mapsto \sin\bigl(\phi_0 + {\textstyle\int_0^t(\omega_0+s(u))du} \bigr) $$
which is the same as
$$ t \mapsto \sin\bigl(\omega_0t + \phi_0 + {\textstyle\int_0^t s(u)du} \bigr) $$

We see that phase modulation is the same as frequency modulation by the derivative of the signal -- which looks like a boost of high frequencies and attenuation of low frequencies.

Your proposed
$$ t \mapsto \sin\bigl( (\omega_0+s(t))t + \phi_0 \bigr) $$
won't work -- it would look like phase modulation, with a modulation strength that increases monotonically and without bound depending on how long it is since $t=0$.


How to compute correct phase in FFT even after applying phase unwrap and zeroing round off error? I am converting a time-domain synthetic 1D signal to frequency domain using MATLAB fft. In the frequency domain, the amplitude vs frequency plot is coming reasonable which means it is showing the correct amplitudes at the desired frequencies. But in the phase vs frequency plot, things are not coming any close to being reasonable.  It is showing incorrect phase values at the frequency values present in the fabricated signal. I have pasted my code here:
Ws=100;  % sample @ 100 Hz

T=10;    % collect data long enough for at least a couple cycles

N=T*Ws;  % compute number samples needed for above

t = linspace(0,T,N);  % and generate time vector

y = (2.2*cos((2*pi*4*t)+(pi/3)))+(1.4*cos((2*pi*10*t)+(pi/6)))+(0.8*cos((2*pi*16*t)));  % and a sample signal

subplot(1,3,1)

plot(t,y)

xlabel('Time(s)')

ylabel('Amplitude')

Y = fft(y)/N;                % FFT

PSD = 2*abs(Y(1:N/2+1));     % and the PSD one-sided

f = linspace(0,Ws/2,N/2+1);  % compute freq vector for Ws

subplot(1,3,2)

stem(f,PSD)                   % plot the result

xlabel('Frequency(Hz)')

ylabel('Amplitude')

Y2=Y;   %store the FFT results in another array

%detect noise (very small numbers (eps)) and ignore them

threshold = max(abs(Y))/10; %tolerance threshold

Y2(abs(Y)<threshold) = 0; %mask out values that are below the threshold

phase=rad2deg(atan2(imag(Y2),real(Y2))); %phase information

subplot(1,3,3)

stem(f,phase(1:length(f)))                   % plot the result

xlabel('Frequency(Hz))')

ylabel('Phase')

I am also attaching the plot which I got. Any help would be greatly appreciated.

There's an error in your script. Your sampling frequency is not 100 Hz like you intended but is actually 99.9 Hz.  Therefore you don't actually have a whole number of cycles in your FFT which creates some leakage. This leakage will also impact the phase estimation.
If you want 1000 points at 100 Hz, the last point should not end at 10 seconds but at 10 seconds - 1/100 = 9.99. You should be careful when using linspace to define a time vector.
Use this n = 1 : 1000; t = (n-1)/100;


What does it mean to say that "sinc" filters are ideal? I suggested to someone that one can create a sinc filter that is ideal enough to be indistinguishable from the ideal, given certain limitations with the data in the first place.

Is that true?

It was only a suggestion yet it was met with unhelpfulness.

An ideal filter (system) is the one which can be described mathematically but cannot be implemented (realized) physically.

Coming to a sinc filter, the ideality of this filter results from its frequency domain definition which is an ideal low-pass filter with zero ripple in the pass and stop bands and zero transition width.

Having these specifications, an ideal lowpass filter, can be described mathematically but cannot be realized using any techniques. However you can still find the resulting impulse response of an ideal lowpass filter by using the formal frequency to time transformation of the frequency reponse of the ideal filter.

This kept in mind, therefore, an ideal low-pass filter is seen to have an impulse response in the form of a sinc(x) signal, which extends from minus to plus infinity in time, which is another manifestation of the ideality of the filter being designed.

When such a filter is to be implemented in practice, it can only be approximated. The approximation being better as the length of the filter increases which however is a very undesired condition.

Instead of forcing approximations to ideal sinc filters this way, it is practically more effective to use other techniques such as weighted windowing to realize those filters, which results in the conversion of filter characteristics from those of the ideal one to the realizable one.


Extract same in phase signal from 3 different sources with shared delayed noise? If I have 3 signals coming from 3 microphones with noise $n$ delayed like this:

$S_1(t) = x(t) + n(t+\phi_1)$

$S_2(t) = x(t) + n(t+\phi_2)$

$S_3(t) = x(t) + n(t+\phi_3)$

Their shared feature $x(t)$ is in phase and exactly the same. $x(t)$ is a voice and the noise can be anything different from $x(t)$.

EDIT 1: To be more precise, the different noises also share similar features but with different phases:

$n_1(t) = \sum_{0 \to n} m_1(t+\tau_1) + m_2(t+\tau_2) + ... +m_n(t+\tau_n)$
$n_2(t) = \sum_{0 \to n} m_1(t+\tau'_1) + m_2(t+\tau'_2) + ... +m_n(t+\tau'_n)$

$etc$

I approximately know what $x(t)$ is, a voice call like "Hello Robot".

How can I extract and detect $x(t)$ from these three signals?

Depending on how good your approximation of x(t) is, you might want to try this:

Subtract your approx. of x(t) from all $S$
Cross correlate the signals to get $\phi_2$ and $\phi_3$ (choose $\phi_1$ = 0)
You now should have n(t),
Subtract $n(t)$ from orig $S_1$.
You now have x(t)

edit: Please note, that this will only work if your approximation of x(t) is quite good, and the noise you have really is equal in all mics (and not introduced due to electronics etc.)


Matlab : Exploring the effects of co-efficient quantization b = [0.9355,-3.6135,5.3603,-3.6135,0.9355]
a = [1,-3.7339,5.3562,-3.4929,0.8752]

Use MATLAB to quantize the coefficients in the vectors a and b to 16 bits of precision. First, take each vector and divide by the smallest power of two such that, after division, the magnitude of the largest filter coefficient is no larger than 1. Next, multiply the resulting vectors by 215, and round to the nearest integer. (You can use the round function.) Divide the resulting vectors by 215. The resulting values should be no less than -1 and no larger than 1. Finally, multiply these values by the power of two by which you originally divided.

I do not understand what that smallest power of two is! Can someone help me out with the code?
Thanks!

For each of the vectors "a" and "b", you need to find the element with the largest absolute magnitude. (I'm not sure of the specific Matlab code for this; probably something like max (abs (a)).)

Then, you need to pick the smallest number of the form 2N (N is an integer) that is larger than the value determined in the previous step.

From there, the rest of the instructions should be straightforward.


How to compare PSD methods? I am trying to compare several PSD methods to determine the most suitable method for my application. Does anyone know a good approach to start with? I am relatively new to this field of topic.
I could possibly compare the results of the psd methods, the computational complexity and the resolution of the results.

How to compare PSD methods?

Same way you deal with any type of engineering problem of this sort:

Write down requirements of your application. What are the priorities, sensitivities and minimum quantitative requirements/thresholds.?
Write down the properties, features and parameter space of the different PSD approaches that you have found.
Compare #1 and #2, find the PSD method and parameters that meets your requirements the best and/or has the best set of tradeoffs.


Pre-processing EEG signals Good Afternoon Everybody.

I am developing filters for EEG signals in these days. During the Develop I find some problems:

The filter is FIR (band-pass using a convolution of low-pass filter and
high-pass filter) they are 5 filters in different bands and the order is high and the coefficients are very small. I have coefficients of range from 3.8145e-10 to 0.2279, and the filter orders are:

Alpha:     133
Beta:      117
Theta:     124
Delta:     11 (IIR filter)
Gamma:     125

The lines codes (only of the filters) are these:

%%Filtro for alpha, beta and gamma signals is the same
[n,fo,ao,peso]=firpmord([Fapp2 Fass2],[1 0],[0.005 0.005],2);
ba=firpm(n,fo,ao,peso);
[n1,fo,ao,peso]=firpmord([Fass1 Fapp1],[0 1],[0.005 0.005],2);
b1a=firpm(n1,fo,ao,peso);
ba=conv(ba,b1a);
Na=length(ba);

%%Filtro señal theta
[n,fo,ao,peso]=firpmord([Ftpp2 Ftss2],[1 0],[0.005 0.1],2);
bt=firpm(n,fo,ao,peso);
[n1,fo,ao,peso]=firpmord([Ftss1 Ftpp1],[0 1],[0.1 0.005],2);
b1t=firpm(n1,fo,ao,peso);
bt=conv(bt,b1t);
Nt=length(bt);

%%Filtro señal delta
[n,Wn]=buttord(Fdpp2, Fdss2,0.005,60);
[bd,ad]=butter(n,Wn);
[n1,Wn]=buttord(Fdss1,Fdpp1,0.005,60);
[b1d,a1d]=butter(n1,Wn,'high');
bd=conv(bd,b1d);
ad=conv(ad,a1d);
Nd=length(bd);Nd1=length(ad);

This part is only the pre-processing part but the coefficients are very
small. the Fs is 250Hz and the delay group will be 0.5320 seconds. If I
change IIR filter to FIR the order is near 200 and the delay group will be
major (I tried that using only the band-pass filter without the convolution but the order was major than the actual orders).
my questions are:

Are there methods for reduce order of the filter without lose quality
as the ripple or attenuation?
If I change FIR by IIR, could I lose or distort information?
Is normal that the coefficients be very small or how could I convert
them in bigger numbers (near to order of 0.01 or 0.001 as minimum)?

Also, I tried to use only High-pass filter and Notch Filter then to do the features extractions but I have similar problem with the delay group too because I want to implement in a FPGA or DSP (for real time process) but the delay group is high compare with anothers works like my work.

Sorry, my english is not good.
Thank you so much for your time.

The advantage to using FIR-filters here is that you can (relatively easily) account for the delay.
Assuming filters with constant group delay (see examples here,) - which seems likely given that your filters have an odd number of taps - all you have to do is to remove (n-1)/2 samples (where n is the filter order) from the beginning of your data after it has gone through the filter.

For example:
Filter your data with the Alpha filter.
It is now longer by 133 samples.
Remove (133-1)/2 samples from the start of the filtered data, and the delay is gone.
Do the same for all of your filters, and you can keep the output of all of your filters in sync.

With an IIR, the delay is not that simple to account for.  Getting the outputs of you filter back in sync would be challenging if you used IIR filters.

If keeping the filtered signals in sync isn't important, you could switch to using IIR filters for all of your filters.

This would greatly reduce the computational load as well as the amount of memory needed.

As with the FIR, there are design methods for IIR that use the desired passband, ripple, and attenuation to calculate the needed coefficients so that you should be able to generate coefficients for the IIR filters that match your FIR filters in everything except having a constant delay.

As for the coefficients, I wouldn't worry about the actual values unless you have trouble representing them on your selected hardware.  Use what your design tool gives you for the coefficients.

There are design tools which can help improve your FIR filters.  Check some of the online tools (like this one.)  Given the separation of the various bands (wikipedia entry for EEG,)  it looks to me like your current filters can't be provoding very good performance, any way.  I get only 22dB of attenuation using 133 taps for the alpha filter.

IIR filters that sharp may be difficult to attain, and may be unstable.

It all depends on how much attenuation you need.  Since the bands are defined with only a 1Hz transition, any high degree of attenuation will require high order IIR filters or very long FIR filters.

As with any design project, it comes down to making the best compromise between performance and cost - and you need to know what performance is required.


Time invariance and linearity of recursive system I am confused about the definition of linearity and time invariance of recursive system

Given $$y[n] = y[n-1] - y[n-2] - x[n]$$

To test time invariance, we shift the input $x[n]$ for shifting and see if the output $y[n]$.

Now, does shifting the input include the recursive terms?

That is if we shift the input by k, then we get
$$y[n-1+k] - y[n-2+k] - x[n+k]$$

and when the output is shifted we get

$$y[n+k] = y[n-1+k] - y[n-2+k] - x[n+k]$$

and thus the system is time-invariant?

I have the same question with regard to linearity.

The system is both linear and time-invariant.

Theory says linearity means that for an input $$ax_1[n] + bx_2[n]$$ your output would be $$ay_1[n] + by_2[n]$$ where $y_1[n], y_2[n]$ are the outputs of $x_1[n], x_2[n]$, respectively.

In your case, an input $x_1[n]$ gives output $$y_1[n] = y_1[n-1] - y_1[n-2] - x_1[n]$$ same for an input $x_2[n]$, thus $$y_2[n] = y_2[n-1] - y_2[n-2] - x_2[n]$$

Multiplying the first equation by $a$ and the second by $b$ we get $$ay_1[n] = ay_1[n-1] - ay_1[n-2] - ax_1[n]$$ and $$by_2[n] = by_2[n-1] - by_2[n-2] - bx_2[n]$$

Adding together, we get
$$\begin{align}
ay_1[n] + by_2[n] &= ay_1[n-1] - ay_1[n-2] - ax_1[n] + by_2[n-1] - by_2[n-2] - bx_2[n] \\
ay_1[n] + by_2[n] &= ay_1[n-1] + by_2[n-1] - ay_1[n-2] - by_2[n-2] - ax_1[n] - bx_2[n])
\end{align}$$
which is the given system with $x[n] = ax_1[n] + bx_2[n]$ and $y[n] = ay_1[n] + by_2[n]$. So it's linear.

For time-invariance, we first delay the input by $k$, thus we get $$y[n] = y[n-1] - y[n-2] - x[n-k]$$

Can you show that for a delayed output (with same delay $k$), the same result is obtained?


Efficient frequency offset estimation of known signals Let's say that my receiver will be receiving known data signals that range from the very simple (tone) to the complex (known data sequence that is modulated (PSK/QAM/OFDM/etc.)).  What are the most computationally efficient ways to get a reasonably accurate estimate of the frequency offset?

Time aligned
If the signals are time aligned, you can conjugate-multiply the received signal with the reference signal and look for the dominant frequency.
Max-Likelihood estimation of the frequency shift and complex gain of a known signal in Gaussian Noise
Let $\vec r$ be the noisy observation of a frequency translated version of a sent $\vec x$ with unknown gain.
Here, let us consider the channel $\mathbf A \in \Gamma$
$$
\vec r = \mathbf A \vec x + \vec w  \in \mathbf C^N \\
[\mathbf A]_{nn} = \rho \exp(\omega n) , \rho \in \mathbf C^N , \omega \in \mathbf R \\
[\mathbf A]_{mn} = 0, \forall m \ne n  \\
p(\vec w|\tau) \propto \exp\left( \frac {-\| \vec w \|^2} {\tau}  \right)
\\
\mathbf A^* \triangleq
\operatorname*{argmax}_{\mathbf A}
\max_{\tau}
\underbrace{p(\vec r| \mathbf A,\vec x ,\tau)}_{  \propto \exp \left( \frac {-\| \mathbf A \vec x - \vec r \|^2} {\tau}  \right) }
\\
\mathbf A^* =
\operatorname*{argmax}_{\mathbf A}
\underbrace{
-\| \mathbf A \vec x - \vec r \|^2
}_{
\triangleq f(\mathbf A)
}
\\
f(\mathbf A)
=
- (\mathbf A \vec x - \vec r )^H
(\mathbf A \vec x - \vec r )
\\
= - \vec x^H
\underbrace {\mathbf A^H\mathbf A }_{|\rho|^2 \mathbf I
}
\vec x
+ 2\mathop{re}
(\mathbf A \vec x)^H \vec r
-
\underbrace
{\vec r^H \vec r}_{
\|\vec r\|^2}
$$
Notice the first and third terms have no dependence on frequency $\omega$ nor the phase of $\rho$.  So we simply need to find the frequency that best aligns $\vec x,\vec r$.  That will imply the phase.
$$
\omega^* = \operatorname*{argmax} \sum_n |\exp(-j\omega n) \bar {x_n} r_n|^2
$$
One can use an FFT or some such frequency estimator to recover the dominant frequency
Not Time-Aligned
If the signals are not time-aligned, then you've got a time and frequency search problem on your hand.  With low SNR, this is essentially the problem solved by the cross-ambiguity function (CAF).  If you've got several dB of margin you can take a shortcut.
If you have a reference signal with some AM variation; you can time-align first by AM detecting both signals and then correlating your input by the reference AM.  If your  reference signal doesn't have much amplitude variation, but its spectrum is "spikey" you might be able to correlate the magnitude spectra of the received and reference.  This is the same trick applied in the frequency domain, rather than time domain.
With the exception of the CAF, these are all fairly cheap operations.


Trying to plot frequency response of a filter Transfer Function in MATLAB, it looks wrong I have the transfer function below which is for a IIR filter and trying to plot its frequency response for omega<pi. r=0.99 and theta=pi/3:

Here is my attempt at implementing it in MATLAB:
N=5000;
W=(0:pi/(N-1):pi);
r1=0.99;
theta=pi/3;
z= exp(j*W);

num= (1-r1);
den= (1-2*r1*cos(theta)*z.^(-1)+r1^(2)*z.^(-2));
TF=num./den;
TFdb= 20*log10(TF);
plot(W,TFdb)

Here is what I get:

Apparently it should look like this:

I am very new to this and can't tell what I'm doing wrong, would appreciate if anyone can help.

How did you derive the second figure? I got the different result using the following code:
b = 1 - r1;
a = [1, -2*r1*cos(theta), r1^(2)];
figure; freqz(b, a)

)
The only problme I can find in your code is that
TFdb= 20*log10(TF);
should be
TFdb= 20*log10(abs(TF));
Did you take the wrong sign of $a_1$ when calculating second transfer function?
b = 1 - r1;
a = [1, 2*r1*cos(theta), r1^(2)];

The above code generates the second figure.


what is the significance of pilot carriers in ofdm in an OFDM system having a radix two FFT (Say N=64) . all of these sub-carriers are not utilized. while some of the carriers are nulls to prevent leakage in to adjacent bands, say a four to five carriers will be pilots. and the remaining will be the data carriers. what is the significance of pilot carriers in OFDM system

regards,
phani tej

OFDM pilot symbols are used to measure channel conditions and/or obtain synchronization information.


Phase shifts of QPSK I'm honestly lost on creating timing diagrams for bit sequences. I understand for QPSK there are symbols 00, 01, 10, 11 with phase shifts of 45,135,225 and 315.

I'm just not understanding how you get these diagrams.
I graphed an example on paper but it doesn't seem consistent (using the bit sequence below).

Maybe I'm a bit lost on the whole even and odd bit sequences.

For example, I'm looking at a bit stream of 01100101 0 (also I assume that lone bit at the end can be written as 00)

So b(t) = 011001010

be(t)=X0100

bo(t)=X1011

E denoting Even, O denoting Odd

Thanks

Edit: Also again, for example, I grabbed this from Wikipedia

https://en.wikipedia.org/wiki/Phase-shift_keying#Bit_error_rate_2

I feel I'm just missing something that will make it all click. Thanks!

In QPSK, the bit information is embedded in the phase of the carrier signal. For example, for bits 00, the phase shift is 45degree. This means the carrier signal $x(t) = cos(2\pi f_c t + \phi)  0 \lt t \lt T$, the phase shift $\phi = \pi / 4$. So $x(t) = cos(2\pi f_c t + \pi / 4) = cos(2\pi f_c t) cos(pi/4) - sin(2\pi f_ct)sin(pi/4)$. The Inphase component corresponding to $cos(2\pi f_c t)$ has factor of $cos(pi/4)$ which is what you see in the I component. Similar logic for Q component corresponding to $sin(2\pi f_c t)$. If the bits are 10, then I component will have factor of $cos(3 \pi/ 4)$ which will be $-\frac{1}{\sqrt{2}}$ so you will see phase shift of 180 degree due to the -1 factor.


Index of stationarity of a time domain signal During my internship, I've come across a software that calculates a percentage called "Index of stationarity". This index shows how much stationary is the signal. Does anyone have any idea on how to assess the stationarity of a signal?

This a very complicated question, and I would say a still open topic. The concept of stationarity is manifold, from pure statistics to applied DSP (strict, strong, wide-sense, quasi-stationarity, cyclo-stationarity, to refer to a recently closed question). The lack of access to sound models and faithful realizations renders the quest quite difficult.

Non-stationary is a non-property (like non-linearity, non-causality). Refocusing on some breeds may help its characterization, whether  you are interested in:

randomly stepping piece-wise constant signals
"sines" of  varying frequency (like chirps)
stochastic data of changing distribution
etc.

Yet, a common framework resides in using time-frequency representations, assuming some local stationarity and trying to assess the non-stationarity according to some metric/distance/divergence (norms, Kolmogorov, Kullback, Jensen, etc.) over the time. Here are two possible (first) references:

Stationary index for abrupt changes detection in the time-frequency plane, 1998
Testing Stationarity With Surrogates: A Time-Frequency Approach, 2010


Why does fine CFO cause rotation of the constellation? Just thinking about it, why does fine CFO cause rotation only? I only see diagrams that show rotation of the constellation.

We assume this is a complex magnitude frequency domain plot
In order for rotation to occur, the complex magnitude would have to stay the same when CFO occurs but the complex argument differ, but clearly, if fine CFO occurs here, it will be sampling a point on the envelope that isn't at the magnitude of one of the points on the constellation after gain equalisation, but inbetween. Especiallly as the magnitudes of the subcarriers can differ if using 16QAM like below.
This does not reflect the typical constellation I see:

CFO is modeled like this $x(t) = s(t) e^{j2\pi f_{\text{CFO}}t} $. This is a phase rotation only and does not effect the magnitude. Convince yourself of it by breaking $e^{j2\pi f_{\text{CFO}}t}$ up as it magnitude ($= 1$) and its phase ($=2\pi f_{\text{CFO}}$). This is all happening in the time domain. If you plot the CFO-free signal and the CFO-corrupted signal magnitude over the top of each other, you should see little to no difference.
This is different from plotting the CFO-free spectrum and the CFO-corrupted spectrum. From the Fourier properties, multiplying by a complex exponential in time will result in a shift in frequency, so the CFO-corrupted spectrum will be shifted by $f_{\text{CFO}}$.
Edit
After seeing the comments, I think I understand the confusion better now. I'm using the term "sampling" to refer to taking a continuous time signal and sampling it. In a receiver, you do it at the symbol period to get a nice clean constellation back.
From what I can tell, the OP uses sampling to mean something else (not wrong just something else). In OFDM, the receiver does a FFT, which is also sampling but in frequency. The CFO introduces a shift in frequency so when the receiver goes to do the FFT, it is going to sample at points where you are getting a component from the main subcarrier and also other subcarriers (ICI).
To answer the question, CFO does not effect the magnitude of the signal (only the phase) and it does effect the magnitude of the spectrum at the subcarrier sampling points.


How is a continuous spectrum for the DTFT possible? So we that a complex sinusoid of the form $e^{j\omega_0n}$ is periodic over $N=2\pi/\omega_0$ only if $\omega_0$ is a rational multiple of $\pi$, otherwise the exponential is not periodic. (see EDIT!)

Then we got to the DTFT, specifically the synthesis equation:

$$x[n]=\frac{1}{2\pi}\int_{2\pi}X\left(e^{j\omega}\right)e^{j\omega n}d\omega$$

Here $X\left(e^{j\omega}\right)$ if a function that's continuous over $\omega$, and $X\left(e^{j\omega_0}\right)$ is how much a complex sinusoid $e^{j\omega_0n}$ contributes to $x[n]$.

But that confuses me, because I feel like not all sinusoids should be able to contribute to $x[n]$ since not all complex sinusoids are periodic. Or are aperiodic sinusoids actually the key to the DTFT?

EDIT: I should clarify that in that first bit where I talk about complex sinusoids being periodic, I'm talking about discrete complex sinusoids, $x[n] = e^{j\omega_0n}$ where $n\in\mathbb{Z}$. For $x[n]$ to be considered periodic there must be some $N\in\mathbb{Z} \mid x[n+N]=x[n]$. For such an $N$ to exist, $\omega_0$ must be a rational multiple of $pi$.

My confusion comes from the fact that the synthesis equation includes values of $\omega$ that are not rational multiples of $\pi$, so aperiodic discrete sinusoids are contributing to synthesizing $x[n]$.

The reason the DTFT has values at frequencies that result in discrete aperiodic sinusoids is because we are interested in representing all possible discrete time signals including (linear combinations of) discrete aperiodic sinusoids. Let's say that we have a signal $x[n] = e^{j \omega_0 n}$ where $\omega_0$ is not a rational multiple of $\pi$. Would you like this signal to have a DTFT? Of course. It would not have one if we did not allow the existence of the DTFT at all frequencies (including irrational multiples of $\pi$).

Put another way, periodicity is not a requirement for the existence of the DTFT of a signal.


Filtering odd and even frequencies I consider a signal of length $N = 2^n$ for some $n$. I want to derive two signal from it, one containing only the odd frequencies and one only the even frequencies. Each of these signals have length $N/2$.

I use this code in MATLAB:

n = 3;
N = 2^n;

signal = ifft( rand(1,N) );

for i = 1:length(signal)/2
    sig1(i) = (signal(i) + signal(i+N/2)) ;
    sig2(i) = (signal(i) - signal(i+N/2)) ;
end

stem(0:N-1,  abs(fft(signal)),'k');
hold on
stem(0:2:N-1, abs(fft(sig1)), 'b*');
hold on;
stem(1:2:N-1, abs(fft(sig2)),'rs');

What I should get are two signals, each one of length $N/2$. One containing only the even frequencies and the other one only the odd ones. The way I proceed can be seen either as the first step of the butterfly algorithm in the FFT, or, by using a FIR filter approach, by applying two different filters to my original signal.

However, the signal containing the odd frequencies is always off and when I compare the FFTs they do not match. Why? What is the explanation? I am a bit puzzled by this result.

I think you need to change your code as follows

n = 3;
N = 2^n;
signal = rand(1,N);

fft_sig = fft(signal);
fft_sig1 = fft_sig(1:2:end);  % even/odd frequencies
fft_sig2 = fft_sig(2:2:end);

sig1 = ifft(fft_sig1);        % the actual time-domain signals
sig2 = ifft(fft_sig2);

stem(0:N-1,  abs(fft_sig),'k');
hold on
stem(0:2:N-1, abs(fft_sig1), 'b*');
stem(1:2:N-1, abs(fft_sig2),'rs');


How the increasing of sampling frequency in OFDM didn't cause increase on the required channel bandwidth In OFDM, how the increasing of sampling frequency didn't cause increase on the required channel bandwidth (fs>>BW), as i know sampling frequency means number of samples per second, if this number of samples increased by logic the the needed bandwidth increase.

You can say the null subcarriers have their own bandwidth if you define an alphabet including "zeros" and use null carriers to transmit these "zeros". As the comment of MBaz,

  If we start with the discrete-time OFDM symbol s[n], and increase the
  sampling rate (defined by the time interval between successive samples
  of s[n]), then the OFDM symbol becomes "shorter" and therefore its
  bandwidth will increase.

But it is not the case. The OFDM symbol $T_u = N_{dft} T_s = 1/\Delta f$.

In this bandwidth interpretation, I think the term "virtual subcarrier" is more appropriate.

What you send over by electromagetic wave is data in $N$ active subcarriers $< N_{dft}$, or the frequency positions of null/virtual subcarriers can be used for other systems; because when you revert back the time domain signal to the frequency domain, what you need is the frequency positions of these active subcarriers.

It is easier to think OFDM as FDM. In this case you have $N$ narrowbands of bandwidth $B_{nb}$ and a sample rate $F_s > N B_{nb}$ is required for practical implementations, such as easily-realizable anti-aliasing filters. OFDM is simply one implementation of this system by using $\Delta f = B_{nb}$ and $F_s = N_{dft} \Delta f = N_{dft} B_{nb} > N B_{nb}$.


Hamming code exercise I'm preparing for my exams I stumbled across this exercise. Now I can't seem to figure it out. Anybody can help me ?

The exercise:

For the 8bit data word: 00111001 the check bits are: 0111. Now if the word is read from the memory and the check bits are 1101, what is the dataword?

The text is translated from another language to English.

I started by doing a check :

The entire word is the dataword (8 bits) combined with the check (4 bits).

In this case: 0 0 1 1 0 1 0 0 1 1 1  1 (the check bits on position 1, 2, 4 and 8, with the right bit being the least significant bit)

Then the check itself:

The check consists of searching for all the bit that have the value one from the dataword in the combined word. It can be seen that position 3 7 9 and 10 have the value one and didn't come from the check. These 4 values are written binary and added up. When this is done. the check is found again. (To confirm that this method works)

  3:  0 0 1 1

  7:  0 1 1 1

  9:  1 0 0 1

  10: 1 0 1 0

0 1 1 1

So this is correct (Just to see if this method is correct)

Then I have another check (1101), but how can I find the data word ?

Thx in advance

How I did it was: Use the original data word and fill in the new check bits. :

001111001101

Then again get the positions where the 1s are and convert that number to binary:

1 : 0 0 0 1

3 : 0 0 1 1

4 : 0 1 0 0

7 : 0 1 1 1

8 : 1 0 0 0

9 : 1 0 0 1

10: 1 0 1 0

1 0 1 0

So this means that the 10th bit is wrong

Since I don't have the actual  solution, i do not know if this is correct, bit it fits the way I learned how to do this.


Find the frequency of a pole (Discrete Form) I am studying a Notch filter. I found its Transfer Function (in Z) and its poles. How can I find the frequency of theses poles (which corresponds to the cut-off frequency)?

Thank you very much for your help

Usually a discrete-time notch filter is implemented by an IIR (infinite impulse response) filter, which has a numerator as well as a denominator polynomial. The notch frequency - the frequency where the input signal should be completely suppressed - is determined by the numerator polynomial, i.e. by the zeros of the filter's transfer function, not (directly) by its poles. Of course, a practical design will also place the poles close to the notch-frequency, but the actual notch is implemented by a zero on the unit circle of the complex plane, i.e. by the numerator polynomial.

The numerator of a second order notch filter looks like this:

$$N(z)=z^2-2\cos(\omega_0)z+1\tag{1}$$

You can verify that this polynomial has two complex conjugate zeros at $z=e^{j\omega_0}$ and at $z=e^{-j\omega_0}$, so the notch frequency is given by the normalized angular frequency $\omega_0$.


How to calculate the gain in a bivariate fft in R? In Statistica gain is defined as follows:

  Gain. The gain value is computed by dividing the cross-amplitude value
  by the spectrum density estimates for one of the two series in the
  analysis. Consequently, two gain values are computed, which can be
  interpreted as the standard least squares regression coefficients for
  the respective frequencies.

However, spec.pgram (the engine behind spectrum) in R does not return the cross-amplitude value (as far as I can tell).  How can I calculate gain for these signals?

Example:

makewave <- function(freq,phase,amp,Nsamples=length*samplerate,samplerate,time=Nsamples/samplerate,as.time.series=TRUE) {
  time <- Nsamples/samplerate
  phase <- phase*(2*pi)/180
  wavetimes <- seq(0+phase,time*freq*pi*2+phase,length.out=Nsamples)
  #plot(1:samples/samplerate,amp*sin(wavetimes),type="l",xlab="Time")
  if (as.time.series) {res <- ts(amp*sin(wavetimes),deltat=1/samplerate)} else {res <- amp*sin(wavetimes)}
  return(res)
}
signal1 <- makewave(15,30,1,180,60)
signal2 <- makewave(15.01,60,1,180,60)
signal.union <- ts.union(signal1,signal2)
sp <- spectrum(signal.union,span=c(3),taper=0)
plot(sp)
plot(sp,plot.type="phase")
plot(sp,plot.type="coh")

I am/was shaky on the correct terminology. I have two signals, and I want to compare them to see how much the 2nd matches the first (coh), how much the 2nd lags behind the 1st (phase)

Normalized cross-correlations are good at that.  I don't know what the relevant command would be in Statistica, but in Matlab you would use the command xcorr.

  and the extent to which the 2nd is under/overshooting the peak amplitudes of the 1st.

To answer that in a useful way it would help a great deal to see the signals that you are talking about.  Are you just interested in the peaks, or the signal in general?  If you are talking about the signal in general then you could time-align the signals using the information from your cross-correlation, and then divide one of the signals by the other.  If you are just interested in the peaks you could do the same and then just extract the results at the peak locations.


Underwater Acoustic Positioning I have to develop an underwater positioning system which has to determine the position of a ROV.

Four buoys will be placed on each corner of a swimming pool. Each one of these buoys will be equipped with an hydrophone and a transducer (which emits acoustic signals). Both devices will be connected to a Raspi for performing DSP.

The ROV will also contain the same equipment (transducer, hydrophone, Raspi).

I was wondering which could be the best way to determine the distances between the ROV and the buoys for triangulating the position.

At first glance I have tought with the next ways to achieve it:

Signal Strength: Every buoy will emit a beacon signal with the same power, but at  different frequency for each one. The ROV could determine its position be means of comparing the power received.
Ping: The ROV emits a ping and inicialices a counter, when the buoys receive the ping they respond immediately. When the ROV receives the ping back stops the counter and calculates the distance.

I think that the second approach can be much more accurate. The first method that comes in my mind in order to achieve it is to use a certain frequency for emitting the ping (let's say 2khz). The buoys would be in recepcion mode, awaiting for receiving this ping. They would detect the ping by means of frequency analysis (Goertzel's algorithm). The same algorithm would be used in to ROV for receiving the ping.

The main drawback is that the system needs to be really precise. An error of 1ms means 1.4m of deviation. I don't really know if that is doable and if is the best way to achieve it.

Which method would you use?

I have worked with this kind of a positioning system before (NOT underwater, but still sound based), and I hope I can help. Your first option is probably not the best idea, even in a controlled environment there are many factors which can change your signal's amplitude and therefore making it unreliable information.

I think your second option is a better way to go. It could be even better if you can embed some time information to your transmitted signal as well (a timestamp like GPS). With a timestamp containing information about when the signal was fired, you can calculate the time difference without ping-back, which would save you some distance error. Or you would have two methods to determine the position and you can perhaps use some sort of filtering depending on a set of initial results (to decide which one works better).

Another method that I would like to suggest is an array configuration. If you can make sure that your transmitters will transmit at ABSOLUTELY at the same time (with respect to each other), you can use the phase differences of the received signals in your boat to determine the range and the respective angles(with some ambiguity). So again you wouldn't need to ping back your signal. But in this configuration the absolute positions of your transmitters must remain the same, and the positioning of your receivers in your boat must remain the same as well.

When we used the ping method, our biggest problem was echo. I don't know exactly how echo works underwater but if it exists, I would suggest you to be careful about it. In our project we had only used the ping without any extras that I suggested to you, and I must admit it wasn't very reliable.

Good luck with your project, I hope you do better than us, and I would like to hear about the results.


Where do I start in learning how to demodulate a BPSK31 signal in software? I'm a software developer and a Amateur Radio operator. Everything I know about audio (which isn't saying much) is self taught.

I've learned how to generate audio samples and build wav files, send those samples to hardware output, etc. I even figured out how to generate BPSK31.

Where I'm stuck now is demodulation. I've read around enough to realize that this is a much harder side of the problem than modulation, and I understand intuitively why this is.

My question is: Where do I even start to learn how to do this? I'm not even sure what the right terms to search for in Google are. So many of the results are hardware solutions or new theoretical approaches to the problem that I've gotten lost in the noise. What should I be looking for to figure this out?

This system uses phase shifts to encode characters (it is a differential encoding scheme). You need to design a system that can detect the phase changes, and that is able to synchronize, that is, find the start of a character. Your receiver must also be able to recover from ocassional errors in the phase change detection. Your question is where to start to learn how to do this; my answer is, read the book linked at the end of this page: http://sethares.engr.wisc.edu/telebreak.html. This book doesn't cover this specific modulation but you can learn everything you need from it.


Vibration Sensor - Balancing Machine Im' working on a balancing machine. this machine has two vibration sensors. (pizzo) my question is about that why this machine use two Sensors? this machine has an rotary encoder to detect disk angle.

i trying to read these signal by ADC card and using FFT bandpass filtering to calculate unbalanced mass on disk. for this algorithm i just use one of Vibration sensors. should i use second sensor?

Of course, the machine could be provided with three, four or more sensors, if the goal of the control unit is perform a balancing of the rotation. Actually three sensors should be the proper selection if the rotation is unconstrained in a biaxial plane, and two sensors, if constrained in a uniaxial direction.


choosing an input signal for system identificaiton I know we can use a

-step
-impulse
-sinus
-multisine
-random phase multisine

To identify a system.

But what are the advantages/disadvantags of each type of input?

conceptually, any driving signal with bandwidth that covers the maximum expected bandwidth of the system being identified is sufficient.  you drive input and synchronously measure input and output, use the FFT to compute the spectra of both input and output, divide the output spectrum by the input spectrum, inverse FFT the resulting spectrum and you have the impulse response.  at least if nothing bad happened, like the system went non-linear because the driving signal exceeded some bounds of saturation.  or if the amplitude of the driving signal was so low that the result of the division has measurement and numerical errors.

so the whole idea of a good driving function is one that does not exceed limits of linearity in the device, yet delivers a lot of energy to reduce measurement and numerical error.  a measure of efficacy for a driving signal is crest factor, which is the ratio of the maximum signal amplitude to the r.m.s. amplitude.  the lower the crest factor, the better.  it means more energy driving the device under test while remaining within maximum limits.

the impulse driving function is natural in that what is measured at the output is the impulse response.  making the nascent impulse thin enough to have wide bandwidth might mean that the input energy is small.  you would need to measure the impulse response several times by driving the system with a train of impulses, spaced apart sufficiently (longer than any expected impulse response) and synchronously average the results.  but the impulse train will be zero most of the time, and if you limit the impulse height to prevent saturation, then the driving signal has low power and your wait is longer in collecting data.

a single sine wave is not sufficiently broadbanded.  it will tell you only how that system behaves at the frequency of the sine.  a collection of added sinusoids at various frequencies will give you a better broadbanded result, but might add up to a spiky driving signal if all of the phases are the same.  scrambling the phases can help, but finding the best "random" phases can be ad hoc.

a swept-frequency sinusoid is a popular driving signal and there is some theory behind it, if the swept frequency varies linearly with time.  if log-frequency is swept linearly with time, the preceding theory no longer applies.  the crest factor is $\sqrt{2}$.

another possible driving signal is a bipolar maximum-length sequence, which is $(-1)^{a[n]}$ where $a[n]$ is a binary {0,1} maximum-length PN sequence.  with that signal the crest factor is as low as it can get (which is 1), but there can be unpredictable behavior if there are sufficient non-linearity in the system being tested.


Why Sparse Priors Like Total Variation Opts to Concentrate Derivatives at a Small Number of Pixels? When performing image deconvolution (deblurring), people often make use of priors to get rid of the illness of the problem. One very common prior is total variation, a sparse prior. Compared to Gaussian prior, it is supposed to be better as it concentrate derivatives at a small number of pixels, which fulfill the assumption people have with clear images. So my question is: how can we understand that sparse prior such as total variation concentrate derivatives at a small number of pixels while the Gaussian prior cannot. Thanks.

I will divide my answer into 3 sections.

The Distribution of the Derivative of Images

Take a real world image, any image.
Apply the derivative operator on it (Namely apply the kernel $ \left[ 1, -1 \right] $ on it.
Display the histogram of the filtered image.

I took this image:

The histogram I got is this:

This distribution is very similar to Laplace Distribution.
Clearly, this distribution is "Sparse" namely most of the values are 0 (Or close) and very few are different.
In real world we assume most of the values which are small yet not zero are actually due to the noise.

Total Variation Optimization Problem

Look at the Optimization Function:

$$ \hat{x} = \arg \min_{x} {\| Hx -y \|}^{2} + \lambda {\| G x \|}_{1} $$

Where $ H $ is the Blur Operator and $ G $ is the Derivative Operator.

Now, you can look at it in the Sparse sense or you can look at it as the MAP solution given a Laplace prior for the Gradient.

As written above, it fits very well to real world images.

Sparse Optimization

Look at the Optimization Function:

$$ \hat{x} = \arg \min_{x} {\| Hx -y \|}^{2} + \lambda {\| G x \|}_{0} $$

Where $ H $ is the Blur Operator and $ G $ is the Derivative Operator.

This optimization problem clearly promote Sparse solution (With respect to the Derivative).
Yet, this is a very hard problem to solve.
Hence it was shown that under some circumstances the solution of the problem with $ {\ell}_{1} $ coincide with this solution.

Moreover, by looking at the Unit Sphere of different norms (And the Pseudo Norm $ {\ell}_{o} $, or more correctly Cardinality Function) you can see why the lower the norm the Sparse Solution it promotes.

All in all you have many point of view on the same problem.


Generate impulse from random frequency response This code supposed to generate random frequency response which could be used as impulse...But it fails to pack correctly when going from fft domain to real signal time domain, the other half is left inside the real signal...

clear all,
close all,

mn = 37.2700     %minimum amplitude value in dB
mx = 140.3640    %maximum amplitude value in dB
N = 44100        %length of the impulse

ampr = mn + (mx+mx)*rand(N,1) %random values for amplitude

a=db2mag(ampr)        %convert dB to magnitude
a1=a(2:end)
a2 = flipud(a(2:end))          %create a symetric mirror for negative frequencies
a3=[a1;  a2]            %combine positive and negative frequency amplitude arrays

pmn =-180        %minimum phase
pmx = 180        %maximum phase

phaser = pmn + (pmx+pmx)*rand(N,1)   %random phase values

ph=wrapToPi(deg2rad(phaser)) %convert phase degrees to radians values between pi and -pi
ph1=ph(2:end)
ph2 = flipud(-ph(2:end))    %flip and create mirror phase values for negative frequencies
ph3=[ph1;  ph2]                %combine two mirror symetric halves of positive frequency and negative frequency phase arrays

mk2=real(ifft(pol2cart(ph3,a3)))  %convert polar values to cartesian and finally to real time domain signal
fs=44100                          %assign sampling rate
% mk2 = mk2(1:floor(length(mk2)/2)) %remove the other half caused by negative frequencies
mk2=mk2/max(max(abs(mk2)))       %normalize the signal

audiowrite('RandomGeneratedImpulse.wav',mk2,fs,'BitsPerSample',24)

Your code looks good and functional, the only thing you need to do is to add DC values for both negative and positive frequencies...Make sure you skip the DC component for negative frequencies when working with odd number of samples...

clear all,
close all,

mn = 37.2700     %minimum amplitude value in dB
mx = 140.3640    %maximum amplitude value in dB
N = 44100        %length of the impulse

ampr = mn + (mx+mx)*rand(N,1) %random values for amplitude

a=db2mag(ampr)        %convert dB to magnitude
a2=flipud(a)          %create a symmetric mirror for negative frequencies
a3=[a; a2]            %combine positive and negative frequency amplitude arrays

pmn =-180        %minimum phase
pmx = 180        %maximum phase

phaser = pmn + (pmx+pmx)*rand(N,1)   %random phase values

ph=wrapToPi(deg2rad(phaser)) %convert phase degrees to radians values between pi and -pi
ph2=flipud(ph)*-1            %flip and create mirror phase values for negative frequencies

ph3=[ph; ph2]                %combine two mirror symmetric halves of positive frequency and negative frequency phase arrays

mk2=real(ifft(pol2cart(ph3,a3)))  %convert polar values to cartesian and finally to real time domain signal
fs=44100                          %assign sampling rate
mk2 = mk2(1:floor(length(mk2)/2)) %remove the other half caused by negative frequencies
mk2=mk2/max(max(abs(mk2)))       %normalize the signal

audiowrite('RandomGeneratedImpulse.wav',mk2,fs,'BitsPerSample',24)


How to sync audio and video signals in MATLAB? I have a video clip that has a frame rate of 25 fps (MJPEG compression) and an audio clip with a sample rate of 48kHz. Both of them play for the same duration - 27 seconds.

Is there any standard way to sync the two signals of different sampling rates if I wish to write them to file?

If this is not doable cleanly in MATLAB, I am happy with a some software tool on Linux or MacOS that can do this syncing.

E.g. the VideoFileWriter object can take both an audio sample and a video frame at a single step, however I have about 1917 audio samples corresponding to a single video frame.

Some details:

I have read in the audio signal using [signal, fs] = audioread(filename).
signal is a 1336321 x 1 array.

All the 697 frames in the video are saved as a cell array of RGB images.

Thanks to helpful comments from @MBaz, I managed to come up with a solution:

we can associate multiple audio samples with a single frame using the VideoFileWriter object. This fact and use-case is missing in the documentation.

First, some stats about the audio and video files. The stereo audio samples are in a 2xN array signal. The video frames are in a cell-array of RGB images, frames.

%% Write audio and video to file
%   Write both audio and video samples into a single video file.
%   Multiple audio samples are matched with one video frame.
%   MATLAB problem: Compression is not possible when audio is included.

% Display A/V stats
fprintf('\nVIDEO STATS:');
fprintf('\nNum. of frames: %d', length(frames));
fprintf('\nFrame rate: %.02f fps', annot.video.frame_rate);
fprintf('\n\nAUDIO STATS:');
fprintf('\nNum. of samples: %d', size(signal,1));
fprintf('\nSampling rate: %.02f Hz\n', infoAudio.SampleRate);

Here's the portion that matches multiple audio samples to a single video frame.
Note that there are 697 video frames and 1336321 audio samples. Some basic fiddling to ensure proper matching. The audio subsampling effect is imperceptible in the output AVI files that I tried.

numAudio = size(signal,1);
numRep = floor(numAudio/length(frames));
numDiff = numAudio - numRep*length(frames); % mismatch

if numDiff
    % if length(frames) does not evenly divide nAudioSamples, then
    % subsample audio to match numRep*length(frames)
    selector = round(linspace(1, numAudio, numRep*length(frames)));
    subSignal = signal(selector, :);
end
assert(numRep*length(frames) == size(subSignal,1));

Finally, we use the VideoFileWriter object available in the Computer Vision System Toolbox of MATLAB to write audio and video to file.

videoFWriter = vision.VideoFileWriter(fullfile(shotPath, 'avclip', ...
                                        [num2str(shotNum) '_av_clip.avi']), ...
                                      'AudioInputPort', true, ...
                                      'FrameRate',  annot.video.frame_rate);
for i = 1:length(frames)
   fprintf('Frame: %d/%d\n', i, length(frames));
   step(videoFWriter, frames{i}, subSignal(numRep*(i-1)+1:numRep*i,:));
end
release(videoFWriter);


Upsample - filter - downsample Given a discrete signal $ r(nT_c)$ specified at a rate $T_c$, assume we want to resample with a 3/2 rate change to match the sampling rate of some other signal.

We would like to upsample it by 3, then filter it by filter of length $K$ and then downsample by 2.

Also assume we have the following relationship exists $$ \frac{3}{2}T_s = T_c$$

Step 1: Upsample then we have

$r'(n\frac{Ts}{2}) = \left\{ \begin{array}{ll}
        r(n\frac{Ts}{3}) & \mbox{if $n = 0,3,6$};\\
        0 & \mbox{if OW}.\end{array} \right.$

Step 2: Filter

$\hat{r}(n\frac{Ts}{2})= \sum_{k=0}^{K-1}r(\frac{(n-k)T_s}{2})h_{filter}(k)\,\,\,\,\,\,\,n=0,1....$

Step 3:

$r_o(nT_s)= \hat{r}(n\frac{Ts}{2}-\underbrace{\frac{K-1}{2}\frac{T_s}{2}}{???}) \,\,\,\,\,n=0,1,$

I don't understand two of the above steps

1) why do we need use a filter to resample and

2)  why isnt it the part underbrace $$\hat{r}(n\frac{Ts}{2}) \,\,\,\,\,n=0,1,....$$

Thank you very much.

Your goal is to create a signal with 3 times more samples. The first step is to insert 2 zeros between each sample.
Technically, you now have a signal with 3 time more samples, as you require. However most of the samples are zeros. What you need to do next is to interpolate the 2 samples between each 2 nonzero samples. And this is exactly what the low-pass filter does (or, at least, this one way to look at it).
Digital filters introduce delay to the signal by half the filter's length. I believe this is the reason for the time shifting


What is the advantage of using a Stepped Frequency waveform over FMCW? What is the advantage of using a Stepped Frequency Waveform instead of a traditional linear FMCW ?

If you're defining the stepped-frequency waveform as a pulse train where each pulse has a different frequency, then it indeed does have advantages. One of the main ones being bandwidth requirements, which we'll go over.
Both FMCW and stepped-frequency techniques use frequency-domain range finding via the DFT. They map beat frequencies resulting from the mixing process to a range. For a given LFM chirp of bandwidth $\beta$ and pulse width $\tau$ the mapping is given by
$$f_b = \frac{-2R\beta}{c\tau}$$
The range resolution capability of any waveform with captured bandwidth $\beta$ is
$$\Delta R = \frac{c}{2\beta}$$
In order to get fine resolution, $\beta$ must be large. In a FMCW system, the entire instantaneous bandwidth $\beta$ must be supported for a single chirp. Once you begin exceeding 100 MHz or so, this can become a big problem.
Stepped-frequency waveforms go around this issue by achieving the desired bandwidth by collecting multiple pulses, with each pulse's carrier frequency increasing by some $\Delta F$ (usually linearly). By doing this, the receiver has to only support the bandwidth of one individual pulse of length $\tau$, which is approximately $1/\tau$.
The desired range resolution can be achieved since the bandwidth necessary is now captured over $M$ pulses and is given by
$$\Delta R = \frac{c}{2\beta} = \frac{c}{2M(\Delta F)}$$
As you may already be thinking, there are disadvantages to this. The timing complexity increases and you must spend more time collecting pulses in order to achieve the desired range resolution.
Edit: Follow-up Questions

Does a step-frequency waveform ease the load on transmitter because its easy to generate?

Traditionally yes, it is more difficult to generate high chirp-rate LFM pulses when compared to simple single-frequency modulated pulses since analog devices were used. Back in the day, having oscillators that could provide the required chirp-rate were expensive to design and build. Today were we usually use a DAC, that limitation has been eliminated. You simply have to meet the sampling criteria, handle the sample throughput, some filtering and you're off to the races.

I would like to understand what kind of problems arise in LFM when the receiver has to handle a large BW?

An example of one of the many issues involves the antenna. Antenna's have defined beam patterns for a particular frequency and are designed with this in mind. When using LFM, the frequency changes intra-pulse and so does the wavelength. This means that the antenna pattern is changing as the pulse is being transmitted. For small bandwidths this can be tolerated. Once you start having hundreds of MHz of bandwidth, the antenna pattern can change drastically, which is usually undesirable.

In a stepped-frequency waveform, is the frequency constant within a chirp? If yes, then you would not get the beat frequency by mixing the transmitted and reflected signal?

In the waveform I described, there are no chirps. It is a train of simple modulated pulses, where each pulse has a higher (or lower if downchirping) frequency than the previous pulse. An example looks as follows:

We collect each individual pulse after it is reflected from a target and "stack" them to construct a pseudo-LFM sequence. Due to the mixing process, this will yield a beat frequency just like in FMCW, and we can take the DFT to process range.
There is a trade-off in the $\Delta F$ chosen and the number of pulses. You ideally want  $\Delta F$ to be small so that you approach an ideal LFM pulse, but then that requires more pulses to capture the bandwidth $\beta$.


Designing appropriate filter by choosing values for $A$ and $B$ Let $$y[n] - ABy[n-1] +A^2y[n-2] = x[n] - Bx[n-1] +x[n-2]$$
Where $A$ and $B$ are constants. Also there are two sounds which have $400$ hertz noise and the sampling rate is $8000$ samples/sec. The task is determining $A$ and $B$ in order to eliminate the noise. I've written a MATLAB program which takes $A$ and $B$ as inputs and then plots the frequency response of the filter as output(using freqz command). Letting $A = 0.9$ and changing $B$ from $1.4$ to $1.9$ leads to shifting the band stop region to the left:

Since we want to remove $400$ hertz noise choosing $B = 1.9$ is good because $0.1 \times \frac{8000}{2} = 400$. The main problem is finding appropriate $A$. Let $B = 1.4$ and increase $A$ from $1.1$. Part of magnitude plot will move downward:

Conversely, if $B = 1.9$ increasing $A$ from $0.1$ to $0.9$ causes plot moves upward. If $A = 1$ then magnitude plot is a straight line

It seems there is some sort of usual trade-off situation here. More elimination of noise is achieved if we select larger values for $A$ but this causes more attenuation in other frequencies as well. Is there any mathematical solution for finding the best values for $A$ and $B$? Also how we can design better filters for removing $400$ hertz noise? I mean alternative solutions instead of mentioned IIR filter.

Bottom Line:
$$A <1$$
$$B =2\cos(\omega_n)$$
Where $\omega_n$ is the normalized angular frequency of the desired notch location (in this case for the OP with a sampling rate of 8KHz and notch at 400 Hz this would be $\omega_n = 2\pi400/8000 = \pi/10$), resulting in $B \approx 1.902$ and $A$ is the frequency notch bandwidth parameter; the closer $A$ is to 1 the tighter the bandwidth of the notch.
Details:
The straight-forward approach to eliminating narrow band noise at 400 Hz would be a classic 2nd order notch filter. This is detailed at this post:
Transfer function of second order notch filter
This approach is done by simply placing a zero at the frequency of interest, and by adding a pole close to the zero, but inside the unit circle of course to be stable, we can adjust the bandwidth of the notch by the proximity of this pole to that zero, meaning $|p|<1$, but close to 1.
Using real coefficients, meaning complex conjugate zeros and poles as detailed in the linked post, this results in the transfer function:
$$ H(z) = \frac{1+a}{2}\frac{z^2-2z\cos(\omega_n)+1}{(z^2-2az\cos(\omega_n)+a^2)} $$
Which as a difference equation is exactly in the form that the OP has used and therefore given this is the same filter, answers the OP's question. This is is derived from the z-transform as follows (we can ignore the gain scaling parameter $(1+a)/2$ which for most cases is close to 1, and express in decreasing powers of z (and assuming I didn't make simple algebra errors!):
$$ \frac{Y(z)}{X(z)} = \frac{1-2z^{-1}\cos(\omega_n)+z^{-2}}{(1-2az^{-1}\cos(\omega_n)+a^2z^{-2})} $$
$$ Y(z) (1-2az^{-1}\cos(\omega_n)+a^2z^{-2})= X(z)(1-2z^{-1}\cos(\omega_n)+z^{-2})$$
$$Y(z)- Y(z)2a\cos(\omega_n)z^{-1}+Y(z)a^2z^{-2} = X(z)-X(z)2\cos(\omega_n) z^{-1} +X(z)z^{-2} $$
And then with the inverse Z transform we get the desired coefficients that the OP wanted:
$$y[n] - 2a\cos(\omega_n)y[n-1] + a^2y[n-2] = x[n]-2\cos(\omega_n)x[n-1] + x[n-2] $$
With comparison to OP's expression:
$$y[n] - ABy[n-1] +A^2y[n-2] = x[n] - Bx[n-1] +x[n-2]$$
Results in the following:
$$A = a$$
$$B =2\cos(\omega_n)$$
Where $\omega_n$ is the normalized angular frequency of the desired notch (for the OP's sampling rate and notch frequency this would be $\omega_n = 2\pi400/8000 = \pi/10$) and $a$ is the notch bandwidth parameter ($a<1$, the closer a is to 1 the tighter the notch, refer to the linked post for complete details).


CODECS used for VoIP calls I am trying to self-study audio codecs and I am trying to understand the difference between CODECS that come under AMR-WB, SuperWB and EVS vs. AAC-LD (MPEG-4). When is the former used and when is the latter used. Can one be used instead of the other? Also where does LPC (Linear Predictive Coding) get used? Any guidance would help

Let me put in an unconnected form things about audio codecs.

Audio can be music or speech (discluding ultrasound, sonar etc).
Music is wideband and requires high-fidelity. Speech has a lower bandwidth and does not require fidelity but intelligibility.

A Codec can be lossy or lossless, the choice depends on the source type, the purpose of application and the resources available.

Codecs are also divided into two as waveform coders (such as AC-3, AAC, MP3, OGG vorbis, (lossless) FLAC, ADPCM) and source coders (such as vocoder, LPC, CELP etc).

A waveform coder, compresses (transforms) the data trying to preserve the waveform shape. On the other hand, a source coder tries to compute parameters of a model that's assumed to produce the given data.

Waveform coders are preferred for high quality multimedia audio, broadcast audio and music compression. Source coders are/can be preferred for very low bitrate speech compression for communication systems such as GSM, military radio etc.

As technology advances, bandwidth limitations have disappeared for certain (but not all) applications. In addition, the compression efficiency of lossy waveform coders improved to the point that medium quality speech compression comparable to that of source coders can be achieved.

So, technically speaking,  your VoIP (voice over internet protocol) application can select from a wide range of speech codecs available for different applications needs and available resources. However you must consult to your VoIP protocol standards to see whether it leaves you free to choose your speech codec or dictates a particular type for use. In short, it's more of a standards issue than a technical one.


Recognizing data clustered into shapes I am working on a project in Python to detect and classify some bird song, and I have found myself in a position where I need to convert a wave file into frequency vs. time data. This hasn't been too much of a problem, but to be able to classify the different syllables into groups, I need to write something that will detect when the data clusters into a certain shape. To give you an idea of what the data looks like, here is an image of how the data looks when plotted:

I need some way to get each individual syllable (each shape with a separation on either side) and save them either to a variable or to their own files so that I can run Pearson correlation between them using SciPy.

Also, I prefer Python, but I am open to coding in other languages if you have another way to do it.

Thanks!

Two questions:

1/ Near 8s, we can observe a stable pitch for 100ms or so, then a sudden increase dropping until 8.5s. Does this whole sequence (8s to 8.5s) form a single entity, or do you consider the two stages (stable then decrease) to be two entities?

2/ Do you want to work with or without supervision. Do you know in advance the "patterns" to look for?

If you want to work without supervision (say you have gathered recordings and aim at extracting a "structured representation" from it), your problem is akin, in a first step, to voice activity detection. Just use signal intensity, maybe in conjunction with a "pitchiness" metric (say the ratio of the maximum of the autocorrelation in the brid range, 1kHz - 5kHz here) to detect segments where there is an active strong pitched tone. Median-filter the resulting sequence to smooth it, and then threshold it to get the different segments. Once you have broken down your signal into segments, you can do interesting things with them. For example, you could extract for each of them a pitch trajectory (a sequence with the strongest frequency peak for each FFT frame, or something more robust extracted with a true pitch estimator), use DTW to compute a matrix of pairwise distances between each block, and use a clustering algorithm (k-means, agglomerative clustering) to identify groups of similar pitch patterns (the 8:8.5 and 10:10.5 segments). It is likely that an unsupervised approach will over-segment - for example 7.6:8.5 and 9.6:10.5 will be recognized as the repetition of the two same blocks, while to you they could fundamentally be one single pattern, but you could use something like Sequitur to have one level of structure higher.
If you have a predefined dictionary of "patterns" you want to label your signal with, you'd better follow the kind of approaches used for speech recognition, with the only major difference that speech recognition does not take pitch into account, while in your case pitch is the sole information to be considered! A speech recognition system tackles both the segmentation and recognition tasks in one single FST decoding operation.


Local Histogram Processing clarification In Gonzales' book Digital Image Processing in the section Local Histogram Processing, he writes, for a 2D image:

  At each location, the histogram of the points in the neighborhood is
  computed and either a histogram equalization or histogram
  specification function is computed and either a histogram equalization
  or histogram specification transformation function is obtained. This
  function is then used to map the intensity of the pixel centered in
  the neighborhood. The center of the neighborhood region is then moved
  to an adjacent pixel location and the procedure is repeated.

He then writes:

  Because only one row or column of the neighborhood changes during a
  pixel-to-pixel translation of the neighborhood, updating the histogram
  obtained in the previous location with the new data introduced at each
  motion step is possible. This approach has the obvious advantages over
  repeatedly computing the histogram of all pixels in the neighborhood
  each time the region is moved over one pixel location.

This last statement I am not clear about. Given a neighborhood, say a 3 by 3 patch centered
at a pixel location $(x,y)$, if we shift the neighborhood over by one pixel to the right say, so it is now centered at $(x+1,y)$, in this case the letmost 3 by 1 column from the previous neighbourhood will not be in the new patch, and a new 3 by 1 column on the right of the previous patch will be present in the new patch.

What does the author mean when he says that this approach has the obvious advantages over repeatedly computing the histogram of all pixels in the neighborhood region each time the region is moved one pixel location?

If you move the 3x3 window to the right by 1 pixel, a new column gets introduced to the right of the previous window and the leftmost column gets deleted.

If the histogram component for the k$^{th}$ intensity level is $p_{k}=n_k/n$ where n is the total number of pixels in the neighborhood and $n_k$ is the number of pixels with intensity k.

For new window, $p'_{k} = (n_k+n_{r_k}-n_{l_k})/n$
where $n_{r_k}$ is number of occurrence of the $k^{th}$ intensity level pixels in the added right column and $n_{l_k}$ is number of occurrence of the $k^{th}$ intensity level pixels in the deleted left column.

$$p'_{k} = n_k/n+(n_{r_k}-n_{l_k})/n$$
       $$p'_{k} = p_k+(n_{r_k}-n_{l_k})/n$$


Does the sound volume affects the MP3 file size? Let's say I have an X music track in some lossless format (wav?).

If I convert it to MP3 format with its original volume and get the file "A.mp3".

Then I edit the lossless source and reduce the volume by 50% (in some way without affecting its quality) and convert it to MP3 and get the file "B.mp3".

Will be any size difference between "A.mp3" and "B.mp3"?

Does the volume reduction cause a loss of information?

Depends on the details, but typically not. For the sake of discussion, we assume CD quality for lossless, i.e. 16 bit, 44.1 kHz.

Reducing volume may or may not reduce information. It depends on what number format it's done in. If you convert 16-bit input to, to for example, 32 floating point, no information would be lost. It may still be lost later, depending on what happens to the data down the road
If you do it directly to the 16-bit input data, you have to
basically throw away the LSB (least significant bit), and your
bit depth has been reduced to 15 bit.
This has little or no effect on MP3 encoding. Most encoders start with a target bit rate and then minimize the encoding artifacts within these constraints. If you encode 10s wave file at 320kb/s, the resulting file will be very 3.2 Mbit long. The thing that varies with content is the sound quality and not the file length.
The type of "information" that encoders use are "masking" and "signal to noise ratio". These are largely independent from the overall scale, so shifting volume doesn't affect the encoders much.
Most encoders these days come in two variants: constant bit rate and variable bit rate. With "variable bit rate" you may see the file length fluctuate a bit, but generally not by much. With constant bit rate, the file lengths will always be constant.


Nyquist Stability Test Currently I'm learning about Nyquist Stability Test and I'm few days stuck on one thing and I don't understand that so I wanted look here for help.

Given is transfer function $K\frac{2s-11}{s(s^2+s+1}$. We take that $K=1$.

Line represens the Nyquist contour for K=1.

Than in my book it's written that
that function is stabil when

$\triangle arg(1+L(Iw))=(max(grad(z_L),grad(n_L)-N_{-}(n_L)+N_{+}(n_L)  )\pi$

than they see from this picture above
$\triangle arg(1+L(Iw))=-\pi$

and

$(max(grad(z_L),grad(n_L)-N_{-}(n_L)+N_{+}(n_L)  )\pi=(max(1,3)-2+0)\pi=\pi$

Not stabil.

What I don't undestand is - $\triangle arg(1+L(Iw))=-\pi$? What this could be? What angle? How do I see that from picture?

The open-loop transfer function

$$G(s)=\frac{2s-11}{s(s^2+s+1)}\tag{1}$$

given in your question has the following Nyquist plot:

which is different from the one in your question.

Note that in any case it's important to know how the trace is closed, i.e., how the points $\omega=0^-$ and $\omega=0^+$ are connected. In this case they are connected by a semi-circle with infinite radius in the left half-plane (the dashed curve in the figure below):

[generated with WolframAlpha].
This is a consequence of the pole of $G(s)$ at $s=0$.

The Nyquist criterion is actually simpler than you made it look in your question. The number of clockwise encirclements ($N$) of the point $-1+j0$ equals the number of right half-plane zeros of $1+G(s)$ ($Z$) minus the number of right half-plane poles ($P$):

$$N=Z-P\tag{2}$$

Since $G(s)$ has no poles in the right half-plane ($P=0$), the number of clockwise encirclements equals the number of zeros in the right half-plane. From the Nyquist plot we see that $N=1$, so there is one zero in the right half-plane, which corresponds to one pole of the closed-loop transfer function in the right half-plane. Consequently, the system is not stable.

This is also easily verified by directly computing the poles of the closed-loop transfer function

$$Q(s)=\frac{G(s)}{1+G(s)}=\frac{2s-11}{s^3+s^2+3s-11}\tag{3}$$

roots([1,1,3,-11])
ans =

  -1.2833 + 2.3182i
  -1.2833 - 2.3182i
   1.5667 + 0.0000i

which shows that there is one pole with a positive real part, as predicted from the Nyquist plot.


Extremely low frequency filters I have a data acquisition board(A/D+Digital Signal Processor) and I want to check if a digital high-pass filter(implemented in DSP) at an extremely low cut-off frequency(0.05Hz) is actually working.

If this was a frequency I could generate with a signal generator it'd be easy to check, but 0.05Hz is too low and I can't generate it. How do engineers check this kind of filters?

if this is a first-order HPF DC blocking filter, this is the recommendation i have for it in fixed-point, and i think it would be a good idea in floating point, too.

while (n<something)
  {
  y[n]   =  x[n]  -  w[n]*( 2^(1-N) );
  w[n+1] =  w[n]  +  y[n]*( (2^(N-1))*(1-p) );
  n = n+1;
  }

where p is the pole (in the z-plane and inside the unit circle) and $0<1-$p$\ll 1$.  N is a big integer.  with fixed-point you would perform the w[n]*(2^(1-N)) multiplication with arithmetic shift right (ASR) or the (w[n]>>(N-1)) operation.  in float you can just multiply it and that simply changes the exponent in the floating point value.  the constant coefficient (2^(N-1))*(1-p) should be precomputed and N chosen so that this is between $\frac12$ and $1$.

the -3 dB corner frequency $\omega_0 = 2\pi\frac{f_0}{f_\mathrm{s}}$ (as normalized angular frequency) is a function of the pole p and is:

$$\begin{align}
\cos(\omega_0) &= 1 - \tfrac{(1-p)^2}{2p} \\
 &= 2 - \tfrac12 \left(p + \tfrac1p \right) \\
 &= 2 - \cosh(\ln p)
\end{align}$$

the pole p is

$$\begin{align}
p \ &= \ e^{\operatorname{arccosh}(2-\cos(\omega_0))} \\
\\
 &= \ 2 - \cos(\omega_0) - \sqrt{\big(2 - \cos(\omega_0)\big)^2 - 1} \\
\\
 &= \ 1 - 2\sin\left(\tfrac{\omega_0}{2}\right) \left(\sqrt{1 + \sin^2\left(\tfrac{\omega_0}{2}\right)} - \sin\left(\tfrac{\omega_0}{2}\right) \right) \\
\end{align}$$

so
$$ 0 \ < \ 1-p \ = \ 2\sin\left(\tfrac{\omega_0}{2}\right) \left(\sqrt{1 + \sin^2\left(\tfrac{\omega_0}{2}\right)} - \sin\left(\tfrac{\omega_0}{2}\right)\right) \ \ll \ 1 $$

and the shift or exponent change N is an integer and should be about

$$ N \approx 1-\log_2(1-p) \approx 1 - \log_2(\omega_0) $$

someone can check my math.

this is an improvement by Tim Wescott to my DSP guru trick.  it deals with both the cosine problem and the limitations of fixed point arithmetic.

latest thinking:  i am not sure, but this might work fine for floating-point operation x, y, and w are all single-precision floats (double-precision floats should be no prob).  (1-p) is a coefficient calculated in advance and is very close to zero.

while (n<something)
  {
  y[n]   =  x[n]  -  w[n]*(1-p);
  w[n+1] =  w[n]  +  y[n];
  n = n+1;
  }


Why is convolution or FFT breaking in this case? It is well-known that convolution in real space corresponds to multiplication in the reciprocal space. Yet, I obtain different results using both methods in the case that I'm going to present. I am sure that I made a mistake, but I really can't see where.

The Python code:

%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt

# Create a signal: a gaussian
sig = np.exp(-(np.linspace(-10, 10, 1024))**2)
# Calculate its FFT transform
F_SIG = np.fft.fft(sig)

# Create an FIR filter from its transform
F_KER = np.empty_like(sig, dtype=np.complex128)
F_KER[:len(F_KER)//2] = np.exp(-1j*np.pi/3)
F_KER[len(F_KER)//2:] = np.exp(1j*np.pi/3)
F_KER[[0, len(F_KER)//2]] = np.cos(np.pi/3)
# Get its inverse FFT transform
kern = np.fft.ifft(F_KER) # is actually pure real value
print("Max real value:\t", np.abs(kern.real).max())
print("Max imaginary value:\t", np.abs(kern.imag).max())
print("Therefore 'kern' can be considered pure real.")
kern = kern.real

# From math books,
# multiplication in Fourier/Frequency space <-> convolution in real space
multFFT_sig = np.fft.ifft(F_SIG*F_KER).real # Again, pure real, can be checked
conv_sig = np.convolve(sig, kern, "full")

# Plot both of them
plt.figure()
plt.xlim([0,1024])
plt.title("Convolution vs. FFT multiply")
plt.plot(multFFT_sig)
plt.plot(conv_sig)

# Padding problem? Circular convolution problem?
plt.figure()
plt.xlim([0,1024])
plt.title("Original signals")
plt.plot(sig)
plt.plot(kern)
plt.legend(["Orig. signal", "Conv. kernel"])

plt.figure()
plt.xlim([0,2048])
plt.title("Zero-padded signals")
p_sig = np.append([0]*512, np.append(sig, [0]*512))
p_kern = np.append(kern[:512], np.append([0]*1024, kern[512:]))
plt.plot(p_sig)
plt.plot(p_kern)
plt.legend(["pad. signal", "pad. conv. kernel"])

# Calculate convolution padded terms
plt.figure()
plt.xlim([0,2048])
plt.title("Padded convolution vs. FFT multiply")
p_F_SIG = np.fft.fft(p_sig)
p_F_KER = np.fft.fft(p_kern)
p_multFFT_sig = np.fft.ifft(p_F_SIG*p_F_KER).real # Again, pure real, can be checked
p_conv_sig = np.convolve(p_sig, p_kern, "full")
plt.plot(p_multFFT_sig)
plt.plot(p_conv_sig)

The output:

Max real value:  0.551327165647
Max imaginary value:     2.55232033425e-15
Therefore 'kern' can be considered pure real.

I also tested with scipy.signal.fftconvolve and compared with the analytical answer: the FFT multiplication gives the right answer (clearly asymmetrical peak in blue) while the convolution (in real space or using FFT) is wrong.

Why is convolution not working as intended?

Yes, zero-padding is the issue, however not as you would except...

@Soren is only partially right: you can add zeros at the beginning, middle or end of the signal just as you did, but absolutely not as you wish.

What's the general rule?

  Look at how the signal repeats with and without zeros and choose the
  place that distorts the signal the least.

What does it mean?

In the case of a real, discrete and aperiodic time-space signal, there is a causal relationship between the value at n and the value at n+1. So all these values should be considered as a bloc, adding zeros in the middle would distort it. Here is an obvious reason from your code, you had a proper Gaussian, not something like this:

In the case of your convolution kernel, it's a bit more complicated. What do the last taps (second half) mean?

In the aperiodic DTFT (i.e. linear convolution, as opposed to circular), these mean that the convolved values depends of both [n, n+1, n+2, ...] (as expected) but also [len - 1, len - 2, len - 3, ...] (where len corresponds to the end of the array without padding). This "classical" linear convolution can be simulated by padding the end/beginning of the signal.
In the periodic DFT (circular convolution), the taps of the second half actually correspond to the values [n-1, n-2, n-3, ...] (not the values by the end of the signal).

np.convolve and bunch try to do a "classical" linear convolution while you seem to want a non-classical one (zeros added in the middle). Here is a comparison:

As I understood, the actual kernel you want is the one described in "your kernel".

You can easily get the right kernel by fftshifting:

kern = np.fft.fftshift(kern)

This has the effect of centering your DC (using linear phase shifting) and then you can just pad zeros at the beginning/end as much as you wish and use the standard convolution/filtering tools.

In the end, this was an interesting question: people usually just pad blindly using convolution functions or options from the FFT transform without ever wondering what they are doing. Some even put fftshifts out of habit every time they deal with anything resembling a kernel.


Are all IIR filters unstable in nature? I know this question has been asked here numerous times. But I am still unclear here :(

My book says,

  A linear time-invariant system is stable if its impulse response is absolutely summable.  (G. Proakis)

So that should mean that, for a stable system, the impulse response h[n] goes to zero as n approaches infinity, right? But that filter is FIR filter. So FIR filters are always stable.

So, for unstable filters, the impulse response is not absolutely summable. In another way, the impulse response never approaches zero. Again, for IIR filter, h[h] continues to go on with n i.e. never goes to zero. So, IIR filters are supposed to be unstable.

But from Z-transform, IIR filters have poles and FIR filters don't, correct? And if the poles are inside unit circle, i.e. if their value is below one, then that filter is stable.

So, what's happening here? If IIR filters have impulse response which is not absolutely summable, how come they are stable when their poles are within unit circle? Their impulse response goes on like forever, right? Their output gets feedback. So why they can be stable depending on pole's location or value?

Thanks in advance :)

Note that for stable IIR filters, the impulse response does approach zero as $n$ goes to infinity. It just never becomes exactly zero. However, the sum of the absolute values is finite. Just as an example, take the exponential impulse response

$$h[n]=a^nu[n],\qquad |a|<1\tag{1}$$

where $u[n]$ is the unit step function. The sum

$$\sum_{n=-\infty}^{\infty}|h[n]|=\sum_{n=0}^{\infty}|a|^n=\frac{1}{1-|a|}\tag{2}$$

is obviously finite (for $|a|<1$).

I agree with Hilmar's answer, but I would like to add that everything he said about poles is only valid for causal filters. In the general case it is not necessary for a stable filter to have all its poles inside the unit circle. For stability it is necessary and sufficient that the unit circle is inside the region of convergence (ROC) of the filter's $\mathcal{Z}$-transform (transfer function). So a stable filter can have poles inside as well as outside the unit circle, but in that case it can't be causal. So it's only the additional requirement of causality that forces the poles of a stable IIR filter inside the unit circle.

A similar thing holds for FIR filters. Causal FIR filters indeed have all their poles at the origin, but non-causal FIR filters also have poles at infinity. A simple example is this non-causal moving average filter

$$H(z)=\frac15\left(z^{-2}+z^{-1}+1+z+z^2\right)\tag{3}$$

which has poles at $z=0$ as well as poles at $z=\infty$, but being an FIR filter, it is obviously stable.