Recover Filter Coefficients from Filtered Noise

It's an interesting problem.  What you have there is what's known as a blind deconvolution problem.  These are well known "hard" problems, but not necessarily impossible.  Finding an algorithm to solve it relies on using some prior knowledge you have about the filter or the noise source driving it.

It's an ill-posed problem mathematically, so if there is a solution, there's no guarantee it's unique.  In other words you might get a set of FIR coefficients that seem to fit the data, but it's not guaranteed that they are the "right" coefficients that model the filter well in the real world.

Having said that, I've tried using the EM algorithm and got some promising results.  It's not perfect, but I think it's possible to get a solution.

I've changed the model slightly so that $$S = N*F + v$$where $v$ is Gaussian noise with some unknown variance.  If we find a good solution, this variance should end up being small compared to the signal.

In the EM algorithm you need some hidden variables: in this case, $N$, and some unknown parameters.  Those are $F$, $\sigma^2$ the unknown noise variance, and $p$ the probability of $N_i\in{0,1}$ being equal to one.  This is the prior distribution on $N$.

The trick to it I found was to use the alternative description of the EM algorithm where the hidden variables are described by an approximating distribution $q$.  We want $q$ to approximate the posterior distribution given the observations and the parameters: $p(N|S,F,p,\sigma^2)$.  This almost certainly doesn't factorize as a product of probabilities $q(N_0)q(N_1)...q(N_n)$ so this is an approximation.

In this case $q$ consists of $q_i$ for $i=0..n$ where $q_i$ is the probability that $N_i=1$ under this approximating distribution.

The EM algorithm is an iterative one, where each iteration involves alternately updating $q$ based on the current value of the parameters, then updating the parameters based on the current value of $q$.

I haven't said anything about $m$ yet: it's unknown, but if we make it too big, the extra FIR coefficients should be set to near zero.  These might be at the start or the end, it just depends where the algorithm decides to converge to.

I simulated some data of my own with the same pulse shape as in your synthetic example, with these results.  The true value of $m$ was 200, so I ran the EM algorithm with $m=400$ to see how it did.  I made $S$ 10000 samples long, to be sure I had enough data.  The result came out at about half the amplitude of the original, and as explained above there are some extra coefficients at the start and the end, but shifting and scaling to compensate you can see that it got the shape of the impulse response fairly accurately:

Zooming in on the middle section:

From your file synthetic_data.txt I got something similar:

As it stands it seems to be overestimating $p$ i.e. putting in too many impulses and underestimating the magnitude of the FIR coefficients to compensate.

Now on test_data.txt, the algorithm gave this impulse response:

It's harder to tell how well the algorithm has done this time because we don't know the true impulse response.  One thing we can do it try to reconstruct $S$ from our estimate $\hat F$ and some estimate of $N$.

I've used $q_i$ to estimate $N_i$ by thresholding: if $q_i>0.5$ I put $\hat N_i=1$ as my estimate of $N_i$, otherwise I use $N_i = 0$.

Convolving that sequence $\hat N$ with $\hat F$ gave a decent match to $S$.  I've zoomed in on a typical short subsequence so you can see the detail:

One caveat I have is that the rate of impulses seem a bit high: about 10% of the reconstructed $\hat N_i$ samples were equal to $1$, which is more than your 1 in 20 to 40 estimate.

But overall I think that's doing well for a blind deconvolution problem, so I'd say it's doable with an EM algorithm or something similar.

There are quite a few other algorithms for blind deconvolution, see the references, and one of the others might be better for you.  Variational Bayesian EM would be worth trying, and there are others which might have some advantages in other ways.  It always depends on how much prior knowledge you have, and on other constraints like runtime.

The code below works in Octave: I haven't tried it in Matlab (because I don't have a licence) but it should work in Matlab too.  Very much prototype code, don't rely on it for anything critical without some more work (sanity checking, catching potential divide-by-zeros etc.)  Runs in about 30 minutes on my two-year-old laptop.

Thanks for uploading the example datasets, which was very useful.

function [Fhat,q] = blindDeconvImpulseNoise(filename)

% Load data from file
S = load(filename);
% Renormalize
scaleFac = sqrt(mean(S.^2));
S = S/scaleFac;

% Estimate (assume filter has less than 400 taps)
tic();
[Fhat,q] = estimate(S,400);
toc();

Fhat = Fhat * scaleFac;

end

% Run the estimation: the EM algorithm
% (see ftp://ftp.cs.toronto.edu/pub/radford/emk.pdf)
function [Fhat,q] = estimate(S,m)

L = length(S);

% Initial values
q                      = ones(L+m-1,1)/2;

% Fhat is the estimated filter.
Fhat = zeros(m,1);
Fhat(max(1,floor(m/2))) = 1;

% Initial values of the scalar parameters: noise variance and
sigmasq_hat = 1;
Rhat = 0.01;

last_sig = inf;

starttime = now;
timeout   = 1800;

while 1

  % Estimate hidden variables
  % Start by converting hidden states into
  Fmtx  = getConvmtx(Fhat,length(q));
  Fmtx  = Fmtx(m:end-m+1,:);

  for i=1:length(q)

    ind = max(1,i - m):min(i+m-1,size(Fmtx,1));
    ind2 = max(1,i - 2*m):min(i+2*m-1,size(Fmtx,2));
    ind2 = ind2(ind2~=i);
    Fi = Fmtx(ind,i);
    resid_less_i = S(ind) - Fmtx(ind,ind2) * q(ind2);
    l = (Fi'*Fi - 2*Fi'*resid_less_i)/(2*sigmasq_hat) ...
        - log(Rhat) + log(1-Rhat);

    q(i) = 1./(1+exp(l));

  end

  qmtx                  = getConvmtx(q,m);
  qmtx                  = qmtx(m:end-m+1,:);

  % Estimate parameters
  Sigma                 = (qmtx'*qmtx) ;
  % Matrix diagonal
  Sigma(1:m+1:end)      = sum(qmtx);
  % Store common term Q'*S
  NS                    = (qmtx'*S);
  % Pseudoinverse estimate of Fhat
  Fhat                  = Sigma\NS;
  % Noise estimate, based on error residual
  sigmasq_hat           = (S'*S - NS'*Fhat)/length(S);
  % Pulse rate estimate
  Rhat                  = mean(q);

  if abs(sigmasq_hat - last_sig) < 1e-10
     % No change on this iteration: stop here
     break
  end
  last_sig = sigmasq_hat;

  if 24*60*60*(now - starttime) > timeout
     warning('Timeout!');
     break
  end
end

end

function Nmtx = getConvmtx(N,m)

Nmtx = zeros(length(N)+m-1,m);

for idx=1:m
    Nmtx(idx:end-m+idx,idx) = N;
end

end


How can I generate a sine wave with time varying frequency that is continuous, in C?

Well lets go lol

@AnonSubmitter85 give to you a nice answer, but let me show my  way to do it in matlab, and this maybe can be very easy to port to C:

First I'm creating 256 samples in 500hz sampled at 44100hz take a look how I accumulate the phase and in the end of the first loop I put the phase between the interval 0 and 2pi...

Nice now lets go to the second loop to create more 256 samples of 1000Hz and to do it continuous I'm using the last phase :-), here is the code:

m_phase=0;

signal=[];

f=500;
fs=44100

phaseInc = 2*pi*f/fs;

for i=1:256
   signal(i) = sin(m_phase);
   m_phase = m_phase + phaseInc;
 end

%place the phaser between the 0 and 2pi range
m_phase = mod(m_phase, 2*pi);

f=1000
phaseInc = 2*pi*f/fs;

for i=257:256*2
   signal(i) = sin(m_phase);
   m_phase = m_phase + phaseInc;
end

%place the phaser between the 0 and 2pi range
m_phase = mod(m_phase, 2*pi);

here is the plot of the above code::


why human voice is band limited while it's time limited?

It's true that band-limited signals can't be time-limited (and vice-versa). However, band-limited and time-limited in that context mean "having finite support", i.e. the respective function is identically zero outside a given interval.

In practice, we almost always deal with functions that are time-limited and band-limited, because we never have infinite time and infinite bandwidth available. However, in that context we use the terms "band-limited" or "time-limited" in an approximative sense. I.e., the contribution of the function outside a given interval becomes negligibly small compared to the noise in that interval. We don't require it to be exactly zero.

Using your example of human speech: if you utter a sentence, no matter how long it is, it will be finished at a certain moment. However, the acoustic signal you've produced will be reflected and reverberated, and it will decay exponentially. Soon the reflections will be inaudible and will be hidden in noise, but in theory they go on forever.

When you say that human speech is band-limited to $4$ kHz, then that's actually not true. We often use that frequency as a cut-off frequency because we know that speech remains intelligible if we filter out frequencies above $4$ kHz. Yet, there is a significant contribution above that frequency, and the harmonics continue indefinitely. But again, their contribution becomes negligible and they disappear in noise.


Why does QAM use two cycles instead of one per byte?

It doesn't matter how long you stay on one constellation point. The constellation diagram has NO ties to time whatsoever. The two base vectors $e_x, e_y$ (x-axis and y-axis) represent your real and imaginary part or "basic functions" your signal is made of. For example, if your signal consists of two functions (in this case, QAM, it consists of a real and an imaginary part including the amplitude)

$$
e_x = A\cos(2\pi f_{carrier}t) \\
e_y = A\sin(2\pi f_{carrier}t) \\
$$

If you stay on one constellation point for an hour, that's fine. The constellation diagram only shows you which "constellation" of amplitude and phase you should choose to transmit a specific "constellation" of bits (001 010 110, etc ..).

The question you ask is about the symbol rate. If the symbol rate is half your carrier frequency, then you will always have two cycles before the next symbol will be sent. If your symbol frequency is a quarter of your carrier frequency, you will have four cycles to the next symbol.

I think there are always at least one and a half cycle shown to make it easier on the eye and for the brain to understand because we can follow the sine easily. This wouldn't be the case if you'd change the phase and amplitude every cycle. That being said, there's also the problem of having enough time for the demodulator to detect the phase shifts. This is why the symbol frequency is often a lot lower than the carrier frequency.


Are delta and double delta features needed when classifying with LSTM

"Very easily" is oversimplification, but overall yes, you do not need to feed deltas into LSTM model.


Looking for algorithm to extract line network from binary images within LabView

The go-to algorithm for extracting edges is to transform the image the an appropriate Hough space. Luckily, the Hough transformation for extracting straight lines is one of the best documented.

It looks like these folks may have a Hough transformation implementation for LabVIEW you could use. If not, it's a pretty well documented algorithm and you could try to roll your own or try to integrate an implementation from another language, eg: MatLab, OpenCV.

With the appropriate Hough transformation, you would be able to extract line segments from your image (converting to a binary image is a great start). These line segments may not connect perfectly, depending on how straight your line is and how clear your image is, and you may need to process the vectors a bit more to get the results you desire.


Ideal sampling - question about the 1/T scaling factor

Yes you are right.

The impulse modulation description of ideal sampling operation introduces an amplitude scale of $1/T$ on the spectrum of the sampled signal, which is a periodic replica of the original bandlimited signal spectrum.

Note that as the sampling period $T$ gets shorter, the scale factor gets larger and the separation between repeating replicas gets larger too. So in the limit you will be left with a single copy of the original spectrum with infinite amplitude scale.

The mathematical definition of energy in the sampled signal is always infinite, regardless of the sampling period (except possibly at the limit of $T = \infty$).

This is an ideal operation. For the practical applications what you are interested in is the finite valued samples of the signal at the sampling instants; the sample values (envelope) will be the same irrespective of the sampling period.


Matching an LPC Magnitude Spectrum to FFT Magnitude Spectrum

I have a similar function to compute the Spectrum from the predicator coefficients, but instead of your code that I do not understand, I have the following :

pOutSpec[spec] = mMeanSquareError / (sumR * sumR + sumI * sumI);

I have something that computes the real and imaginary parts into my buffer :

data[2 * spec] = sumR / (sumR * sumR + sumI * sumI);
data[(2 * spec) + 1] = -sumI / (sumR * sumR + sumI * sumI);

Then I use different functions to display the spectrum according to the chosen format (Real, Magnitude, Square Magnitude or Log).

I also have some problems with scaling, both for the FFT Spectrum and the LPC Spectrum. I have chosen to scale them based on my display area, but as a result, their scale is not exactly the same on the display.


Signal Reconstruction after fourier transform

Suppose you want to perform N points FFT, then the evenly spaced frequency vector is given by

  f= (0:NumUniquePts-1)*Fs/N;

where NumUniquePts = ceil((N+1)/2) is the number of unique points in f, and Fs is the sampling rate.

So if fftx=fft(x,N), then $\alpha_1$ = fftx[15/(Fs/N)+1] and $\alpha_2$ = fftx[40/(Fs/N)+1] (provided 15 and 40 can be divided by Fs/N).

See here for more.


Simple FFT convolution with impulse response wav file doesn't give the right results

Your code performs the circular convolution.
Circ Convolution gives an array equal to max( len(a) , len(b) )
Convolution will return an array of length = len(a) + len(b) - 1 // a, b > 0
If you think about it, adding a reverb makes a sound longer. That's the equation ^^ that describes by how much. What you're hearing is the end of the reverb played over the start.
Keep in mind in a musical setting, often the original sound is mixed in with the reverb sound.

I zero padding your data before FFTing it.

long bufSize = 1024*62;
for(int i=bufSize; i<2*bufSize; ++i) res[i].mRe = res[i].mIm = res_ir[i].mRe = res_ir[i].mIm = 0;
mFFT->runFFT(res, bufSize);
mFFT->runFFT(res_ir, bufSize);

double a=0;
double b=0;
double c=0;
double d=0;
for (int i=0;i<bufSize*2;i++)
{
    a=res[i].re();
    b=res[i].im();
    c=res_ir[i].re();
    d=res_ir[i].im();

    res[i].mRe = (a*c - b*d) * (gain);
    res[i].mIm = 0;
}

That gets rid of the 'wrapping' effect of circular convolution at least.
Also, the imaginary part is zero because the convolution of two real signals is real.


image processing find liquid level (only with images processing, no AI)

Constrain the environment. Best is to make sure the glass is always the same (position and shape), take a picture of the glass empty and use that as reference. Otherwise, you'll need some algorithm to first detect the position of the glass (if the background is constrained, you can use that).
If the liquid is isolated, I'd recommend applying a sobelY filter and something to reduce noise. Take a vertical line out of the graph and identify each step (top glass, water level, bottom glass).


How to find transfer function by state space representation matrices

Use the fact that $\dot x_1 = s x_1$ and take it from there. It appears that your result is correct.

It should be noted that there is a closed form expression to map a state-space representation to a transfer function. However, it includes inverting a matrix with s variables and may hence be more challenging than just solving it as stated above


Sum of sines sound wave using microcontroller

To explain the problem (without solving it):

  However, for some reason, besides seeing the pure 17KHz frequency, I also hear (and see on the oscilloscope) parasitic 5KHz frequency!

Your sampling rate is 80 kHz, so the Nyquist frequency is 40 kHz.  Your generation method produces a distorted sine wave at 17 kHz, which produces harmonics at:

17×2 = 34 kHz
17×3 = 51 kHz
17×4 = 68 kHz
17×5 = 85 kHz
...

But everything above 40 kHz gets aliased, so

51 kHz becomes 40-(51-40) = 29 kHz
68 kHz becomes 40-(68-40) = 12 kHz
85 kHz becomes 40-(85-40) = -5 kHz, which then aliases off 0 Hz and becomes 5 kHz
...

  I've also noticed that for frequencies such as 16KHz (index increased by 80) ,8KHz there are no parasitic frequencies.

Now you've got harmonics at 32, 48, 64, 80, ...  These alias to 32 kHz, 32 kHz, 16 kHz, 0 kHz, etc.  So they aren't as obvious, but they're still there unless your wavetable is lining up perfectly with the samples.

If you look at your output using an FFT analyzer, you should be able to see these aliased harmonics more clearly:

http://www.music.mcgill.ca/~gary/307/week4/wavetables.html


Loosing information in decimation - without aliasing

The whole point of aliasing and the sampling theoreme is that you cannot (generally) know what an aliased signal was, as you cannot represent infinite bandwidth using finite infornation.

If you have extra knowledge of the input signal (eg that DC is not possible, then you might deduce what a string of 1-1-1 was originally.


What is the procedure to obtain a 3D face dataset? Can anyone suggest a good 3D face dataset?

I don't personally know of a good one to vouch for, but on face-rec.org/databases/
at least, isn't the "GavabDB: face database" free?

see:
http://gavab.escet.urjc.es/recursos_en.html

has 3 (apparently free) downloads:
Brief description: GavabDB is a 3D face database. It contains 549 three-dimensional images of facial surfaces. These meshes correspond to 61 different individuals (45 male and 16 female) having 9 images for each person. The total of the individuals are Caucasian and their age is between 18 and 40 years old. Each image is given by a mesh of connected 3D points of the facial surface without texture. The database provides systematic variations with respect to the pose and the facial expression. In particular, the 9 images corresponding to each individual are: 2 frontal views with neutral expression, 2 x-rotated views (±30º, looking up and looking down respectively) with neutral expression, 2 y-rotated views (±90º, left and right profiles respectively) with neutral expression and 3 frontal gesture images (laugh, smile and a random gesture chosen by the user, respectively).

Note:
Those publications that use this signature date must reference the following work: A.B. Moreno y A.Sanchez. GavabDB: A 3D Face Database. Proc. 2nd COST Workshop on Biometrics on the Internet: Fundamentals, Advances and Applications, C. Garcia et al (eds): Proc. 2nd COST Workshop on Biometrics on the Internet: Fundamentals, Advances and Applications, Ed. Univ. Vigo, pp. 77-82, 2004

Readme.rtf (9 kB; .rtf)

CD1_faces1-20.zip (38 MB; .zip)

CD2_faces21-40.zip (39 MB; .zip)

CD3_faces41-61.zip (38 MB; .zip)


Analysing 2500 frequencies using FFT with an input vector of 2048 samples?

The first 1024 output bins of the FFT represent frequencies from 0 Hz to Fs / 2 = 2.5 kHz. So each bin is 2500 / 1024 = 2.44 Hz wide. In other words you have a resolution of 2.44 Hz in the frequency domain.


How does spectral density of a power signal change when it is multiplied by $\cos(2\pi f_ct)$?

$$R_{n}(\tau) = E(n(t)n(t-\tau))$$
$$R_{n_i}(\tau) = \cos(2\pi f_ct)\cos(2\pi f_c(t-\tau))E(n(t)n(t-\tau))=\cos(2\pi f_ct)\cos(2\pi f_c(t-\tau))R_{n}(\tau) $$
Using $\cos(a-b) = \cos a \cos b + \sin a \sin b$, we get
$$R_{n_i}(\tau) = \cos(2\pi f_ct)[ \cos(2\pi f_ct ) \cos(2\pi f_c\tau)+ \sin(2\pi f_ct) \sin(2\pi f_c\tau)]R_{n}(\tau) $$
or
$$R_n(\tau) = \frac{R_{n_i}(\tau)}{\cos(2\pi f_ct)[ \cos(2\pi f_ct ) \cos(2\pi f_c\tau)+ \sin(2\pi f_ct) \sin(2\pi f_c\tau)]}$$


Implementation of goertzel algorithm in C

It's been a while, but if anyone would like a working fixed-point implementation (no floats/doubles) in C you might take a look at: http://www.ti.com/ww/cn/uprogram/share/ppt/c6000/Chapter17.ppt
– Slide 14.

It looks like it will work beautifully and efficiently in even a low-powered micro-controller.
(In which case some of the ‘int’ variables in that TI C code may need to become ‘longs’ to be 32-bit.)


Comparing spectra from hyperspectral images

The goal is to predict a genetic distance of species from the spectra. For instance if the genetic distance of species 1 and species 2 is very height I would expect a different spectra. If the distance is 0 or very low,m the spectra are the same /very similar.

In that case, defining the distance between two spectra based on the spectra makes no sense – you want to define it based on the genetic distance.

So, you want to infer a model that describes genetic distance based on spectra. That's fine – but such a model will be very complex, and it's (as far as I can tell) not possible to find an easy function that will just give you what you want.

Unless you have a large dataset of genetic distances, I'll even go as far as saying that this is scientifically a very dangerous thing to do: your spectra are themselves large sets of data (a lot of points, say 1000). If you have, say, only 100 genetic distances, then you'll always find something that correlates well with the genetic distance within your large number of points, without having anything to do with it: the 101. genetic distance would not be represented by that model at all, i.e. it would only fit to the data set you had when you formed that model.

Basically, what you're describing is a problem that these days would usually be solved using neural networks. Again, this all depends on you having enough data to train your network, or else your network will not generalize to new spectral measurements and be useless outside your training data set.


Fourier transform of a discrete periodic signal

For a discrete periodic signal you normally compute the Fourier series coefficients. They can be obtained from the DFT by dividing the output of the DFT by N, where number of samples in a period.

The DFT assumes the signal is N-sample periodic anyway.

In Matlab you will implement this as

x=[-3 1 1];
ck=fft(x)/3;


How do I compare chords played on a guitar with a known sample?

You probably need more than one algorithm for an accurate detection. Work on frequency domain (mostly via DFT and spectrogram) is the most often used initial transform. After that sophisticated probabilistic models originally used in speech recognition are applied, such as hidden Markov models, dynamic Bayesian networks, and conditional random fields.

The review article Automatic Chord Transcription Using Computational Models of Musical Context gives you a thorough overview on the relevant technique, although it focuses on a more complicated chords recognition with more than one instrument (key, bass note, beat are all considered as well).

There are also some open source chord detection algorithms:
Chordino, Chordata, and
LabROSA Chord Recognition. Note that in LabROSA, a supervised learning method is applied to "teach" the model recognize the chords with higher accuracy.

This is a very interesting topic, and I look forward to see your update in your guitar chord detection result.


You detect the chirp, and subtract it from the mixture. The remainder is your signal of interest (plus other sources of noise which you don't mention). Then you calculate the power of each component, and calculate the ratio. That's SNR.

The challenge is usually detecting the chirp, but that can be done, depending on the level of information you have about it (chirp rate, phase of the chirp, beginning of the chirp, bandwidth, amplitude) using simple techniques like a correlation with a prototype chirp, or "blinder" methods like letting a PLL run on the signal and watching when it locks.


Get histogram of edge lengths/orientations in an image using OpenCV

I don't know if OpenCV has some implementation - but i can provide you with a pointer of algorithm for this.

MPEG 7 provides a fairly well defined (and well thought out) edge histogram descriptor [EHD] which essentially serves your purpose.

You can refer to this: Efficient Use of MPEG-7 Edge Histogram Descriptor which i think is an excellent paper on this.

If you want to see some implementation of MPEG-7 you can find it on JOANNEUM RESEARCH MPEG 7 page.

For more information on MPEG-7: List of possible image features for content based image retrieval


Distinguishing Spoken Speech from Played-back Recording of Speech

Do any general-purpose techniques exist to distinguish live spoken
  speech from played-back speech samples?

No, at least not how I understand your definition.

Consider the case of  "live speech" : some one talks into the microphone, you digitize the signal and analyze it. You can also store the digitized speech and analyze it a day later. Now the exact same samples constitute "play back". Since the samples are identical, you can't distinguish them from "live speech".


remove zero padding effect crosscorrealation

so let's say the length of your FFT is $N$.  let's say that you fill half of the buffer with your signal and fill the other half with zeros.

$$ x[n] = 0 \quad \text{  for  } \tfrac{N}2 \le n < N  $$

the DFT of $x[n]$ is $X[k]$.  then magnitude-square and inverse DFT:

$$ R_x[n] = \mathcal{iDFT}\bigg\{ \bigg|X[k]\bigg|^2 \bigg\} $$

in the time domain, this comes out as

$$ R_x[\ell] = \sum\limits_{n=0}^{\tfrac{N}2 - 1 - \ell} x[n] x[n+\ell] $$

we call $\ell$ the "lag" of the autocorrelation.

now suppose $x[n]$ is periodic with some period $P$ which is much smaller than $\frac{N}2$.

$$ x[n+P] = x[n] \quad \text{for } 0 \le n < n+P < \tfrac{N}2  $$

with a lag of zero

$$ R_x[0] = \sum\limits_{n=0}^{\tfrac{N}2 - 1} (x[n])^2 $$

which, if we apply a little bit of statistical sleight-of-hand, we might expect is $\tfrac{N}2$ times the expected power $\overline{x^2}$.  in the limit with $P \ll \tfrac{N}2 $we can sorta do this with periodic $x[n]$.

$$\begin{align}
R_x[0] &= \sum\limits_{n=0}^{\tfrac{N}2 - 1} (x[n])^2 \\
       &\approx \sum\limits_{n=0}^{\tfrac{N}2 - 1} \overline{x^2} \\
       &= \frac{N}2 \overline{x^2}
\end{align} $$

now with a lag of $mP$ where $m>0$ is an integer:

$$\begin{align}
R_x[mP] &= \sum\limits_{n=0}^{\tfrac{N}2 - 1 - mP} x[n] x[n+mP] \\
       &= \sum\limits_{n=0}^{\tfrac{N}2 - 1 - mP} x[n] x[n] \\
       &= \sum\limits_{n=0}^{\tfrac{N}2 - 1 - mP} (x[n])^2 \\
       &\approx \sum\limits_{n=0}^{\tfrac{N}2-1 - mP} \overline{x^2} \\
       &= \left( \frac{N}2 - mP \right) \overline{x^2} \\
\end{align} $$

eliminating $\overline{x^2}$ we have

$$ R_x[mP] \approx \left( 1 - \frac{2mP}N \right) R_x[0]  $$

for $P \ll \tfrac{N}2$.

this suggests an envelope on the autocorrelation of

$$ \big| R_x[\ell] \big| \le \left( 1 - \frac{2 |\ell|}N \right) R_x[0]  $$

this is for the rectangular window.  it can be generalized for other windows.  the envelope results from the two windows having less area in common as the lag $\ell$ increases.


Can anyone help with cross correlation error in Simulink?

The edited code in your question should not result in an all-zero vector. Of course, since x contains many zeros, the autocorrelation also has many zeros, but it's non-zero in its center portion (i.e., around index $2000$).

Also note that you don't need to define a second identical vector to compute the autocorrelation. You can use the command z = xcorr(x);


Palmprint Identification - Why do we align the images before we use the Fourier Transform?

From the article:

  [...]only a sub area of a palmprint image is used in feature extraction.

The alignment by a shift and other image transformations is done in the full hand image which is then cropped to a sub image containing only the features of interest. To analyze the sub image they Fourier transform it. They are not shifting the sub image alone but the full image so the shift operation does not have a direct equivalent operation in the frequency domain of just the sub image.


What is the set of most predictable or in other words most self-correlated signals?

Sums of discrete cosines are good candidates because they're bound to repeat after a finite number of samples. In fact,  anything that repeats after some finite number of samples is the most predictable thing in the world because you can predict the entire thing from just a handful of samples. In fact, the converse of this statement is true, i.e. anything that repeats after a finite number of samples can be written down as a sum of discrete cosines with different phases. In fact, so can your own candidate:
$$c = \sum_{k=0}^{\infty}c_k\cos(\omega_k + \phi_k)$$
where $c_k=c\delta(k)$, $\omega_0=0$, $\phi_0=0$. In other words a constant signal is just a cosine of zero frequency. So there's no tautology here.

You may say of course, that your example is still better, but if we're trying to establish classes of predictability, then both your example and the generic form fall into the same class. I don't see a contradiction here.


How to find the instantanous SNR over Fading channel

You are trying to find the average SNR. In instantaneous SNR, the power of the channel is $|h|^2$.

Your SNR equation isn't accurate. $N_0$ is usually denote the PSD of the noise. The noise power is $\sigma_n^2 = N_0\,W$, where $W$ is the signal bandwidth. Or you define SNR as $E_b/N_0$ where $E_b$ is the signal energy.


First of all, sampling frequency and sampling rate are synonymous. You mean sampling frequency of 44.1 kHz with a data length of 20 ms.

Second of all, 20 ms of data at 44.1 kHz will give you 882 points. Not enough for a 1024-point FFT.

You will either need to upsample to 51.2 kHz, that way 20 ms will give you 1024 points. Or you could append 142 zeroes to your data. Or, you could simply compute your FFT with the 882 points that you have. All 3 methods, will give you the same information. I personnaly recommend that you simply compute your FFT on 882 points.

As to the actual syntax, just google it, there are some pretty good examples.

edit : This seems like a good place to start

https://www.oreilly.com/library/view/elegant-scipy/9781491922927/ch04.html


Why should an image be blurred using a Gaussian Kernel before downsampling?

An image "should not be blurred using a Gaussian Kernel" in general. This can be a safe bet for a lot of basic image processing needs, and a smoothing is almost mandatory when you want to control the information lost by the downsampling.
Blurring is (often, not always) another word for low-pass filtering. When an image contains high-frequency content (fast variations), downsampling without blurring can produce visually weird, annoying aliasing artifacts or moiré patterns. There is an example on wikipedia Aliasing: the original

and the downsampled, represented at the same size:

The ripples on the bottom right are low-frequency artifacts generated by a careless brute-force downsampling. A blurring would attenuate image sharpness, dim the borders between bricks, and reduce the apparent aliasing aspects.
The choice of appropriate blurring filters has a long history in image processing. Gaussian shapes have long been considered somehow optimal for different reasons for "theoretical" continuous images. Plus, it is both decreasing and symmetric in  the space and the frequency domains. In the time domain, this means that faraway pixels have less influence. In the frequency domain, frequencies are reduced monotonously from low to high.
Since most images are discretized, reality is somewhat different. Since the Gaussian convolution  used to be computationally expensive, early approximating filters were designed with short support, borrowed for instance from Pascal triangle. Later, fast recursive implementations were designed (Deriche, Shen, etc.)
I guess question 1) is answered. For question 2) simple averaging gives an equal weight to all pixels in the window. Hence, faraway pixels are given equal importance with respect to closer pixels, which is not optimal in regions where images exhibit weak stationarity, like trends, edges and textures.


What's the computational complexity of MMSE equalizer in OFDM system

For a SISO channel, assuming you have estimated channel coefficients for each sub-carrier $k$, the MMSE equalizer is
$$
\hat{h}_k=\frac{h_k^*}{|h_k|^2 + \frac{\sigma_x^2}{N_0}}
$$
So you can see already there is a multiplication in the numerator ($h_k^*$), and then there is a division by a term in the denominator. So for all $N$ subcarrier, this itself will be $2O(N) \rightarrow O(N)$ for large $N$. There is additional $|h_k|^2$ computation which is basically multiplication with conjugate. So if you have precomputed SNR ($\frac{\sigma_x^2}{N_0}$) and $h_k$ for all sub-carriers, complexity should still be linear on $N$. But this is a very simplified assumption because in most practical OFDM packets, there a pilot sequence periodically present so you may need to update $h_k$ using available pilot sequence.


how to plot fft with imaginary and real part in 3d or calculate the degree of rotation

(1) nparange -> np.arange, and
(2) np.fft.rfft(len(fftSinus),d=1/1000)
You want the frequency values. You can start with just np.arange(1000) for normalized frequency, or frequency indices, and worry about actual values later.
Organize your ax.plot() as follows:

x = frequency index
y = imag
z = real


Plotted ECG signals are not around Amplitude 0 line

I have tried the BEADS technique we have developed for the separation of

sparse signals, possibly positive, with sparse derivatives,
a trend (low frequency),
random noise.

I have played with one signal for about two minutes, and here is one outcome:

It is called BEADS, for  "Baseline Estimation And Denoising with Sparsity". The Matlab BEADS code is at MatlabCentral, the paper is Chromatogram baseline estimation and denoising using sparsity (BEADS). It was initially meant for chromatographic signals, but is also used for EEG/ECG and several other spectral signal, see the BEADS page .

One nice feature is that you directly see the noise (and how uncorrelated it can be) and the shape of the trend, then tweak the parameters. Mine (probably not optimal) are:

addpath('C:\Users\duvall\Documents\MATLAB\toolbox\2014_BEADS_Baseline\');
% load('s0010_rem.mat')
 data = val(1,:)';
% Filter parameters
fc = 0.004;     % fc : cut-off frequency (cycles/sample)
d = 1;          % d : filter order parameter (d = 1 or 2)
% Positivity bias (peaks are mostly symmetric)
r = 1;          % r : asymmetry parameter
% Regularization parameters
amp = 0.8;
lam0 = 0.5*amp;
lam1 = 5*amp;
lam2 = 4*amp;

[x1, f1, cost] = beads(data, d, fc, r, lam0, lam1, lam2);
figure(2);
subplot(4,1,1)
h=plot([data]);axis tight;grid on;set(h,'LineWidth',2);
subplot(4,1,2)
h=plot([x1]);axis tight;grid on;set(h,'LineWidth',2);
subplot(4,1,3)
h=plot([f1]);axis tight;grid on;set(h,'LineWidth',2);
subplot(4,1,4)
h=plot([data-x1-f1]);axis tight;grid on;set(h,'LineWidth',2);

Basically, it took me 10-times more  to download and load the signal than to tune the parameters (90 seconds approximately), but I know the code :)


Refining accelerometer noise using Kalman filter

but what are the steps when only one sensor is available and the physical model of the moving object is not available

Then it's not a Kalman filter.

A Kalman filter works because the system is observable.  In hand-wavy terms, you need to have redundant information about your system states, either because you have actual redundant inputs, or because you have an adequate* model of the system dynamics, and you're watching the system output over time.  If you can't do those two things somehow, then you don't have the prerequisites for a Kalman filter.

Taking GPS/IMU fusion as an example, your design process can be completely naive about the actual plant dynamics because you can model the plant as a rigid lump of stuff that takes rotation rates and accelerations as inputs, and gives geolocation as outputs.  Your model has states for position, velocity, and orientation, and because you're constantly correcting these against the GPS readings, you can get a pretty darned accurate estimate** of the states.

As another example -- and possibly as an inspiration for you -- I'm very slowly developing a Kalman filter for the motion of an object that's constrained to a 2-sphere (in less hoity-toity mathematics, that's the surface of a sphere in 3-dimensional space).  Because we have gravity, it'll know "down"; as soon as the object turns on the sphere it'll know its altitude -- but it'll never know the compass direction***.

So -- if you can somehow come up with redundant information, then you can possibly come up with a true Kalman filter.  This redundancy can be as tenuous as the one I'm using in my second example, where the fact that the object is always a fixed distance from a point tells me a lot.  It could be something like a button on an app that says "ok, Phone, I promise that I'm holding you still" -- I'm pretty sure that with a 6-degree of freedom IMU you could measure horizontal and vertical distances that way (i.e., you could make a tape measure app), but I doubt you could do more.

* What is adequate?  Enough to make the system observable.  Which is a circular definition within the context of this answer, because it would take too much space to explain.  It's a concept that you'll find about two weeks into a senior/graduate class in state-space control systems.  The bottom line for this answer is that it's a property of the system model -- and you don't have one of those.

** Assuming you have a rich enough motion profile -- in order to get good orientation information, the system needs to see varying acceleration, and the corresponding changes in position.  Such a system, if it was bolted to a rock, would be able to estimate position, velocity, and "down" very well, but unless it had gyros with errors significantly smaller than the Earth's rotation speed it wouldn't be able to tell you North and East.

*** Unless I use magnetometers, or get that really good gyro I mention above.


why did convolution coding achieve BER less than Reed-Solomon?

Reed-Solomon codes are good for correcting bursts of errors (such as those that often occur in computer memory systems such as RAM, or on less volatile media such as CDs and DVDs) or at the output of block decoders -- either a received word is completely correct, or it has multiple errors all within the codeword length -- or at the output of convolutional decoders -- if a decoder goes off the correct path, there is a burst of errors until the decoder gets back on the correct path. Reed-Solomon codes are not all that good for correcting random errors such as those that occur in AWGN channels.  It has been common practice for many students (as well as their faculty member advisors, who really ought to know better), not to mention engineers in the communications industry, to put up graphs comparing how well convolutional codes do compared to Reed-Solomon codes in preparation for touting the advantage of their thesis work or product over the competition.
As a matter of fact, systems needing to meet the most stringent performance criteria often use a concatenated coding system in which the outer code is a Reed-Solomon code and the inner code (the one whose codewords are transmitted over the channel) is a convolutional or turbo or LDPC code, all of which are well-matched to the characteristics of the channel.


For filters with impulse responses on a specific $x[n]$, how do you calculate DFT of $y[n]$?

How would you find Yk (the coefficients of the DFT of y(n))?

By applying the DFT to Y.

So: math.

Mathematically, you filter something with a FIR by convolving the signal with the taps vector:

$$y = x * h$$

Hence, the DFT of $y$ becomes is:

$$\mathcal F \{y\} = \mathcal F \left\{ x * h \right\}$$

and thanks to the properties of the DFT, convolution in time domain becomes multiplication in frequency domain:

$$\mathcal F \{y\} = \mathcal F \{ x\}  \mathcal F \{ h \}$$

So you just apply the DFT to both your signal and your taps vector and multiply.

Seriously, this is very basic DFT theory, and you will learn that if you continue your studies of the discrete Fourier transform.


You should not be using floats or doubles to store your message. You need to write your compressed message in binary.


A (Clear) Line Between Computer Vision and Image Processing

I believe Gonzalez and Woods are competent enough to rely on their opinion:

  There is no general agreement among authors regarding where image processing stops and other related areas, such as image analysis and computer vision, start. Sometimes a distinction is made by defining image processing as a
  discipline in which both the input and output of a process are images.We believe this to be a limiting and somewhat artificial boundary. For example, under this definition, even the trivial task of computing the average intensity of an image (which yields a single number) would not be considered an image processing operation. On the other hand, there are fields such as computer vision whose ultimate goal is to use computers to emulate human vision, including learning and being able to make inferences and take actions based on visual inputs. This area itself is a branch of artificial intelligence (AI) whose objective is to emulate human intelligence. The field of AI is in its earliest stages of infancy in terms of development, with progress having been much slower than originally anticipated. The area of image analysis (also called image understanding) is in between image processing and computer vision.

So I would say, that primary difference is in goals, not methods. For example, if the goal is to enhance image for later use by humans, than this may be called image processing. And if the goal is to emulate human vision (be it object recognition, defect detection or automatic driving), then it is closer to computer vision. Note, however, that emulating human vision by definition may also require image enhancement, so in most real cases computer vision relies on image processing.

Image understanding (feature extraction) may be equally used in both - pure image processing and computer vision.


frequency components of discrete Gaussian-windowed signal

No window will give you zero processing (coherent) gain. Check out the table in fred harris's paper.

The table says there is a processing gain of $0.37$ to $0.51$ depending on your $\alpha$ parameter:


Nonuniform sampling points based on the first order derivative of signal amplitude

A method that I use very often to approximate the CDF is based on the observation that given a set of $N$ samples, the CDF at the point x the fraction of the points that are smaller than $x$. So if have a vector of sorted samples $S$ you could approximate the CDF as $cdf(x) = i/N$ if $S_{i-1} < x < S_i$.
A simple way to plot the CDF would be
plot(sort(S), (1:N)/N)
Then you can approximate the derivative by a difference equation. You could use a binary search to find the region of the vector $S$ where $x$ lies, say $idx(S, x)$, then you could compute something like
$i_1 = idx(S, x - \Delta)$, $i_2 = idx(S, s + \Delta)$ and use $(N \cdot (S_{i_2+1} - S_{i_1})) / (i_2 - i_1 + 1)$.
A different way to choose the anchor points would be, let $i = idx(S, x)$ and compute $(N \cdot (S_{i+k} - S_{i-k})) / (2k+1)$
And you can fit use a polynomial fit if you need a better accuracy where the CDF derivatives varies rapidly.


Clarification on how bandwidth effects input signal

The system $H(f)$ appears to be a lowpass filter. As you state, the fact that the bandwidth decreases results in a more distorted output. However, I would not say that "fewer segments of the signal pass through the system".

The Fourier transform of the input in your image is a $\mathrm{sinc}()$ function. This function extends to $\pm \infty$, so to get absolutely no distortion you would need an infinite bandwidth. As this is not possible, then distortion will always appear. Depending on where we "cut" the sinc's "tails", the output will be more or less distorted. If we let "more frequencies" pass, then the output will look more like the input.

The square pulse has very sharp edges - it suddenly goes from $0$ to $1$. This speed translates to high frequency components. As $H(f)$ is a lowpass filter, high components will be filtered out. To what extent this will affect the signal will depend on the bandwidth. As it decreases, more high frequencies will be cut off, leading to a slower transition at the output, which leads to that smoothed version of the square pulse.


This is the most simple FT that you should understand. The basic rule is: a sine of frequency f0 in time domain is a delta at f0 and -f0 in frequency. You can do this theoretically, as shown in http://mathworld.wolfram.com/FourierTransformCosine.html.

Those peaks in your graph are consistent with this. However, it isn't true that there is a peak in at 90Hz. The peaks should be at 10 Hz and -10 Hz. There are two possible approaches here:

a) Use fftshift to correct this issue. This will take each frequency to its real place.

Fs=100; %sampling frequency
T=1; %signal length
N=T*Fs;     %number of samples
f=-Fs/2:Fs/N:Fs/2-Fs/N; %frequency vector
x=cos(2*pi*10*t);
Xf=fft(x);
Xf=fftshift(Xf);
plot(f,abs(Xf)) %magnitude

b) Assume the signal is real, so the spectrum is symmetric. It is true in this case, but it doesn't have to. You could just plot the first half of your fft, i.e. from 0 Hz to 50 Hz.

By the way, if you plot plot(f,abs(Xf)/N), you will find that the peaks have an amplitude of 0.5A (being A the amplitude of the sinusoid). The phase plot seems OK to me, especially after fftshift-ing it.


It's not particularly important since the different bases just result in a different overall scale factor, but the basic shapes are the same.

This being said, any actual numbers need to be defined with a specific bases otherwise it's ambiguous. The most common is the $dB$ which (for power) is defined as
$$L = 10 \cdot log_{10}(P/P_0)$$

i.e. it uses base 10.


What happens to signal, if we offset periodic auto-correlation?

The periodic autocorrelation function is not defined correctly in the OP's question: it should be
$$c_n = \sum_m x_m x_{(m+n) \bmod N} ~~\text{or} ~~ \sum_m x_m x_{(m-n) \bmod N}$$ depending on taste and not $c_n = \sum_m x_mx_{(n-m)\bmod N}$ which is a convolution not a correlation operation.  Note also that there is an
implicit assumption that $\mathbf x = (x_0, x_1, \ldots, x_{N-1})$
is a real-valued sequence, not a complex-valued sequence. Of course,
the autocorrelation sequence $\mathbf c = (c_0, c_1, \ldots, c_{N-1})$ is also a real-valued
sequence with the additional property that $c_0 \geq \max_n |c_n|$

With this correction, consider the sequence
$\mathbf d = (c_0-A, c_1-A, \ldots, c_{N-1}-A)$,
that is, $d_n = c_n-A$ for $0 \leq n < N$.
Let the DFT of $\mathbf c$ be
$\mathbf C = (C_0, C_1, \ldots, C_{N-1})$ where the $C_n = |X_n|^2$ are
nonnegative real numbers.  Note that the DFT of any
autocorrelation sequence cannot have negative values.
Now, the DFT of $\mathbf d$ is obviously
$\mathbf D = \mathbf C - (NA,0,0,\ldots, 0)$ and so we see easily that

  $\mathbf D$ cannot be the DFT of an autocorrelation function
  if $D_0 = C_0 -NA < 0$, that is, if $A > N^{-1}C_0$.

Note that $\displaystyle X_0 = \sum_{n=0}^{N-1} x_n$ and
$\displaystyle C_0 = |X_0|^2 = \left|\sum_{n=0}^{N-1} x_n \right|^2.$
Let $\displaystyle \bar{x} = \frac{1}{N}\sum_{n=0}^{N-1} x_n = \frac{1}{N}X_0$.

  So, let us assume from this point on that $A \leq N^{-1}C_0 = N\left|\bar{x}\right|^2$.
  Is there
  a sequence
  $\mathbf y = (y_0, y_1, \ldots, y_{N-1})$ with
  autocorrelation function $\mathbf d$ whose DFT is $\mathbf D$?

The answer is Yes,  Consider a sequence $\mathbf y$ with
$n$-th term $y_n = x_n - a$. Its autocorrelation function
is
$$\sum_m y_m y_{(m+n)\bmod N} = \sum_m (x_m-a) (x_{(m+n)\bmod N} - a)
= c_n + Na^2 - 2aX_0.$$
We want the quantity on the right to have value $c_n-A$, and so if
we can find an $a$ such that $A = 2aX_0- Na^2$, we are done.
Solving for $a$, we get that
$$a = \frac{2X_0\pm \sqrt{ 4\left|X_0\right|^2-4NA}}{2N}
= \frac{1}{N}X_0 \pm \sqrt{\left|\bar{x}\right|^2-A/N}
= \bar{x} \pm \sqrt{\left|\bar{x}\right|^2-A/N}$$
Note that the quantity in the radical is nonnegative since
we have assumed that $\frac AN \leq \left|\bar{x}\right|^2$
and so we do obtain a real-valued $\mathbf y$ as long as
$A$ does not exceed its maximum value $N\left|\bar{x}\right|^2$.
When $A$ does have its maximum possible
value, $\mathbf y$ whose
$n$-th term is $y_n = x_n - \bar{x}$ is just $\mathbf x$ "centered"
to have "DC value" $0$ since we are subtracting off the average
value of $\mathbf x$ from each component.


With Frequency Shift Keying, the modulation (digital data) takes up bandwidth, so you can't just keep only the frequencies of the mark and space tones. A firm lower bound on how little bandwidth you can use is the distance between the mark and space frequencies, plus half the baud rate on either side. So for 1200 baud with frequencies of 1300 hertz and 2100 hertz, the absolute minimum bandwidth is (1300-(1200/2)) [700 hertz] to (2100+(1200/2)) [2700 hertz] which is a bandwidth of 2 kHz. People have tried to filter it tighter but if the reception still provides the correct data, it is only because of chance. Usually there is also some pulse shaping in the FSK signal before modulation to make the filter's job easier.


The xcorr function assumes a linear cross-correlation, i.e. it assumes that the signals are zero outside the intervals. If you want to have an efficient implementation of the periodic cross-correlation, you can refer to the properties of the Fourier transform with

$$
\mathcal{F}(u \star y)=\mathcal{F}(u) \cdot \mathcal{F}(y)^*,
$$

while exploiting the fact, that the DFT considers signals to be periodic anyway. Here's the corresponding code for that:

N = 20;

u = randn(N,1);
y = randn(N,1);

R = intcor(u, y);

R2 = real(ifft(fft(u) .* conj(fft(y))))/(N);  % Calculate the periodic CC

hold off;
plot([u y R])
hold on;
plot(R2, 'ko');

legend({'u', 'y', 'intcor', 'fft-based'});

To answer your questions in the comments:

It is still not clear to me, if I can use xcorr to calculate the autocorrelation of a periodic signal correctly? -
No, you cannot use the xcorr function, as it assumes a linear cross-correlation. What you want is a calculation, where you assume the signal is periodic, and you just input one period of the signal.
What is the difference between linear cross-correlation and periodic cross-corellation? - Consider the general equation for cross-correlation:
$$
R_uy(h)=\sum_{k=-\infty}^{\infty}u(k)^*y(k-h)
$$

Normally, the summation is done over all times (i.e from $-\infty$ to $\infty$). However, the signals you pass to the xcorr-function are time-limited (i.e. they consist only of $N$ non-zero samples). Now, in the linear case, the cross-correlation assumes that the input-sequences vanish for $k<0$ or $k>N-1$. In the circular-correlation (i.e. assuming periodic signals), the signals are not assumed to be zero for the index above, but instead we assume $u(k+iN)=u(k), \forall i\in\mathbb{Z}$. I.e., the sequences are assumed to be periodic with period N.


python: how to compute the sharpness features of image

What you have implemented is not what the formula denotes, but what you have implemented will return an indication of Sharpness, possibly with some minor modifications.

The formula produces a single $\mathbb{R}$eal number as the ratio of the local sum of laplacians divided by the local mean. "Local" here refers to a square image patch (of some dimensions) centred around pixel $(x,y)$.

If you forget about the division by $\mu$ for a minute, the $LP(x,y)$ is the sum of the "local" second derivatives. This is a $\mathbb{R}$eal number.

The way you have implemented it, the sum operates over the Laplacian of the whole image.

Furthermore, as the Laplacian is a $\mathbb{R}$eal number and the image is a $\mathbb{N}$atural (bounded) number (a.k.a a non-negative integer), you may experience rounding errors. I say "may" because the matrix used by opencv is also based on integers, so the rounding comes into play at a different region.

Evaluating the Laplacian over the whole image is not so much of a "problem" if that sum you used operates across both dimensions. Usually, in platforms such as GNU Octave, MATLAB and others, a single sum over some matrix, operates across one dimension (usually it sums the columns). Therefore, your s would become a vector...not a $\mathbb{R}$eal number.

Hope this helps.


How can create a .wav bar based spectrogram on matlab?

You can get the bar based spectrogram by decomposing the signal segments into fractional-octave subbands. Try octaveFilterBank and then compute signal power level by pow2db(sum(abs(signal))). This will result in a logarithmic frequency scale.
You can also use pwelch to estimate power spectral density, which brings you a linear frequency scale.


Which digital modulation scheme would be suitable for transmitting digital data through the output of a PC soundcard

You've got a pretty good set of circumstances here; you should be able to meet your goal without too much trouble. I don't see anything in your description that would eliminate a whole class of modulation (e.g. phase-shift keying, frequency-shift keying, etc.). Some of the factors that would go into the choice of a suitable format would include:

The required spectral efficiency (i.e. how much data throughput you need relative to the available bandwidth)
The complexity requirements for your receiver (which is usually the most complicated part of the system)
How much effort you're willing to put into developing the implementation.
Other circumstances specific to your application (e.g. if you have poor timing accuracy at one or both ends, known interference, or a poor channel response)

So, ticking these off one by one for your system, we can come up with a few guidelines:

It seems that your biggest constraint is your channel response (which is limited by your sound card's DAC). If you have 40 kHz of one-sided bandwidth available, then you're going to be limited to a symbol rate that is somewhat below that. For a target data rate of at least 40 kilobits per second, you'll want some scheme that transmits multiple bits per symbol.

Provided that your embedded platform isn't loaded down with too many other functions, a modern 120-MHz ARM processor should easily be able to handle demodulation of most any format in the tens-of-kilobits-per-second range.
I'm not sure specifically which model you're running with, but many recent processors provide very tight integration of onboard ADCs with the memory and interrupt subsystem, perhaps allowing you to (without manual CPU intervention) automatically sample the input signal at a specified rate, store the samples in onboard memory, and trigger a processor interrupt only when a certain-sized block of samples is available to be processed. I know that some Atmel devices at a minimum provide this sort of functionality; I've had good success with them in the past.

I assume that this is probably one component of some larger system, so you don't want to have to build something really complex in order to meet your requirements. Luckily, based upon what you've said, I think you should be able to implement something that meets your requirements pretty simply. I assume that since the two devices are cabled together, you're likely to have a very high signal-to-noise ratio (40 dB or greater), which for digital modulations is typically measured as $\frac{E_b}{N_0}$. That makes life much easier, as you likely could skip additions like error-correction encoding altogether (or at least implement a very simple scheme that would catch the very infrequent errors that you would expect at such an SNR).

As far as special circumstances go, for this system, I wouldn't expect much of any. I would expect the oscillator accuracy on the PC side to be pretty good (at a minimum crystal-controlled, so you're in the range of <50 ppm or so; possibly much better if the oscillator is calibrated using some other more precise source). The embedded side is likely to be the same; I assume you're using a crystal oscillator as the clock source. Since the two ends are cabled together, I'll assume that you don't have interference of note.

So, smashing all of these together into a single recommendation, I would probably start down the path of a quadrature phase-shift keying (QPSK) approach at 24 kilosymbols per second. At 2 bits per symbol, this yields a data rate of 48 kilobits per second, which is in excess of your requirement. This particular rate makes your implementation a bit easier; since the output DAC is running at 96 kHz, this results in 4 samples per symbol (it's always easier to run at an integer number of samples per symbol time). I would probably try to design the embedded side so that it would sample at the same 96 kHz rate if possible; this avoids the need to do any resampling on the more resource-starved end.
To avoid any problems with the DC notch that your sound card DAC employs, you could modulate the QPSK signal onto a carrier at 24 kHz. Then, the modulated signal's spectrum would have a null at DC, which would line up with your notch. It's possible that the notch could end up not being a problem at all (especially if it is truly just a few Hz wide as your suggested). In that case, you could potentially get by with an even simpler scheme that just works at baseband, bypassing carrier modulation altogether.
QPSK is a good choice because of its simplicity, both at the transmitter and receiver. At your SNR, you could accomplish greater spectral efficiency using a more complicated scheme like quadrature amplitude modulation (QAM), but the constant-envelope property of PSK signals is attractive from a receiver complexity standpoint. As a note, if you truly needed more bits per symbol in the future, you could move to a higher-order PSK constellation like 8- or 16-PSK. However, these are suboptimal from a bit-error-rate performance standpoint when compared to QAM constellations.
As far as a library implementation goes, I'm not aware of anything that you would just be able to drop in and go, especially for an embedded platform. Your receiver implementation is likely to be tied in with the hardware interface to some extent. You may be able to find some existing implementations for the various steps required for the demodulator, but you'll need to at least tweak what you might find to work well on your platform. The GNU Radio project is a good place to look if you just want to see C++ implementations of many different communications signal processing operations, and it might even yield a useful framework for implementing the transmitter onboard your PC. As a summary, the high-level steps that your receiver would need to do would include:

If a nonzero carrier frequency was used:

Frequency synchronization: locate and track the carrier frequency offset due to oscillator mismatches between the transmitter and receiver (in many cases, that frequency offset will be approximately constant over time)
Carrier demodulation: translate the signal to baseband, yielding two in-phase and quadrature (I/Q) signals (often expressed as a complex baseband signal).

Matched filtering: pass the baseband signal through a filter matched to the pulse shape used at the transmitter (you can likely get away with a rectangular pulse)

Timing synchronization: locate and track the times that correspond to symbol transitions; this can be done by tracking the peak locations in the matched filter output modulo the symbol time

Bit slicing: convert the matched filter outputs at symbol sample times to hard bit decisions

Serialization: make sure you write the multiple bits per symbol out in the right order!

This might sound like a complicated process, but building a practical receiver for even a simple situation like this can be very illuminating. Just comment if there's anything else I left out.


Python: Dynamic Time Warping, what actually is a 'similarity score' and how to compute it?

You have got exactly what you requested in your code.

DTW is a similarity measure between time series. By default, tslearn
uses squared Euclidean distance as the base metric (I am citing the
documentation).

Another ground metric can be used, when specified in the code.
You use a default, i.e., Euclidean-distance metric, your sequences (time series), $s1$ and $s2$, are aligned $len(s1) = len(s2)$. The Euclidean distance is a square root of sum of squared differences of sequence items: $DTW_{Euclid}(s1,s2) = \sqrt{\sum_k(s1[k]-s2[k])^2}$. For identical sequences you have: $DTW_{Euclid}(s1,s1) = \sqrt{\sum_k(s1[k]-s1[k])^2} = 0$. The idea of similarity as the qualitative degree, when referring this Python ecosystem's feature, is supplanted by the concept of distance metric, although the word remains. You have no other choice as to accommodate yourself to this word usage.
This is how tslearn computes the "similarity" when you use a default metric. For other metrics and unaligned sequences, consult the documentation. The portal dsp.SE is ready to assist you in the case of uncertainty.


What is absolute ripple specification and how does it arrive in low pass prototype for Butterworth?

In the first formulation, $\Omega_c$ is the $3$ $\textrm{dB}$ cut-off frequency, because for $\Omega=\Omega_c$ you obtain $|H(\Omega_c)|^2=\frac12$. In the other formulation, the attenuation at the pass band edge $\Omega_p$ is determined by the constant $\epsilon$. For $\Omega=\Omega_p$ you have $|H(\Omega_p)|=\frac{1}{1+\epsilon^2}$. I.e., by choosing the constant $\epsilon$ appropriately, any desired attenuation at the passband edge can be achieved.


How can I objectively measure train noise from a sound recording?

One of the radio stations around here used to have the slogan "It may not be your favorite song, but it has a lot of the same notes."  A DFT analysis can tell you what notes are present, and at what intensity.  With calibration, you can even measure the real world level.  This doesn't get anywhere near psycho-acoustics though.
What you could possibly use it for it to define an objective measureable standard of when the greasing needs to be done.  The loudness of the noise is going to be dependent upon where you measure it, so calibration and absolute values aren't that useful.  What you need, as you alluded to, is comparable recordings of the sound when the tracks are well greased vs when they get annoying.
With those in hand there are several techniques you can use to make quantitative comparisons.  The difficulty for you is that your intended audience will likely not have the technical understanding to comprehend any sophisiticated analysis so you are going to have to try to reduce down to a simple number or two, or a nice visual image that expresses the data succinctly.  I think you already have the latter with your spectogram.  For the former, a simple sum of the magnitudes of the bins in the frequency range of interest should be adequate.
Now, if you can track such a metric throughout a greasing cycle or two, then you can establish objective criteria (by association) of when the noise has reached the same level of annoyance.  Just because you are a concerned citizen, maybe you could email this value to the railroad company for every train passing and they will gladly monitor your numbers, thank you for your service, and send somebody out every time it reaches a threshold value.
Dream on?  Yep.
If there is a local noise ordinance, you can perhaps work with your local politicians  rather than the transport authority.


Impulse response of an LTI system given the input and output signals


Physically understanding autoregressive model (really basic question)

Your understanding is right. Just a few points...

I think you meant Yule-Walker equations.
It is better to distinguish $N$ from $n$. In an AR($N$) model we consider last $N$ terms to approximate the next term.
An important interpretation that you missed was the fact that when we consider
$$\epsilon=x[n]-\tilde{x}[n]=x[n]-\sum_ka_kx[n-k],$$ or in a more standard way (representing the process by variable $y$ instead of $x$
$$y[n]=\sum_ka_ky[n-k]+\epsilon[n],$$
it can be seen as an input-output relationship of an LTI filter, $\epsilon[n]$ being input and $y[n]$ the output. More specifically,

   AR model can be seen as the output of an all-pole IIR filter which is excited by a white noise input.


A little bit better than zeroing bins in order to filter

The Fourier transform of your window-shaped frequency response is its impulse response.  Evaluate the transform of your window function, look at its length above your desired noise floor, adjust the filter response as required to make the impulse length appropriately short, add zero-padding of this length to your original signal before the FFT, and you approach classic convolution filtering.

The added zero-padding is required to remove or reduce artifacts due to circular convolution using FFT/IFFT fast convolution.  Knowing the length of the impulse response is required to know how much zero padding is needed.  Smoother (with wider transition bands) frequency responses tend to produce shorter impulse responses.  And using smooth window shaped functions tends to imply smoother wider transition bands (although there are many other methods that may do better).


how twiddles are calculated in the code snippet below from ti.com?

These are similar and provide the same result. In the c code every other sample:  w[k+1], w[k+3], etc, are the imaginary twiddle factors in the MATLAB code:

For instance the lines in C:

w[k]     = cos (2 * PI * i / n);
w[k + 1] = sin (2 * PI * i / n);

Corresponds to this line in MATLAB :

twiddle(mm+1) = cos(theta) + (1i*(sin(theta)));

Where "theta" in MATLAB is "2*PI*i/n" in c.


Difference of doing a PSD estimate of data and logarithmic transformed data?

The exercise is instructive. Since this is an exercise, I will not do the full problem.  I leave that to you.

grabbing the lynx data from

https://www.encyclopediaofmath.org/index.php/Canadian_lynx_data.

The data is counts and there are no years with a zero count.

and putting it into matlab and taking the log and subtracting the mean of the log,

Would you say that the data looks sine like ?

This technique works on certain kinds of data.  Think of it as a tool and is covered in a number of basic stats book.


What is the $\mathcal{Z}$-transform of a constant?

if you define the constant like this:
$$ x[n] = C ~~~,~~~\text{ for all } n $$

Then its $\mathcal{Z}$ transform
$$ X(z) = \sum_{n=-\infty}^{\infty} x[n]z^{-n} = \sum_{n=-\infty}^{\infty} C z^{-n} $$

does not converge for any value of $z$. Hence its $\mathcal{Z}$ transform does not exist.

Note that $\mathcal{Z}$ transform must be an analytic function with continuous derivatives of all orders unlike the DTFT which is not necesarily an analytic function and which would admit the use of impulses (or non-continuous functions) in its expressions.  Then this would allow the DTFT of the constant $x[n] =C$ to be $X(e^{j\omega}) =C 2\pi \delta(\omega)$. But as stated, impulses (which are not analytic functions) are not allowed to exist in $\mathcal{Z}$ transforms.

Had you defined the constant as the practical signal like
$$ y[n] = C u[n]$$
where $u[n] = \begin{cases} 1 &,& \text{ for } n \geq 0 \\ 0 &,& \text{ for } n < 0\\ \end{cases} $ is the unit step function, then its $\mathcal{Z}$-Transform would exist:
$$ Y(z) = \sum_{n=-\infty}^{\infty} y[n]z^{-n} = \sum_{n=0}^{\infty} C z^{-n} = \frac{C}{ 1- z^{-1} } ~~~, ~~~ ROC~~ |z|>1$$

The distinction is clear I suppose.


What is the definition of half duplex?Is this system,$A$ give $x$ to $B$,then $B$ give $x$ to $C$,called half duplex?

The receiver/transmitter pair in B's possession is operating in half-duplex mode.

As you say "$A$ give $x$ to $B$, then $B$ give $x$ to $C$" which specifies that while $B$ ls listening to what $A$ is telling it (the $x$ that is being given to $B$), the transmitter in $B$'s possession must remain silent until $B$ has completely received $x$. It is only then that $B$ can start up his transmitter and send $x$ to $C$. So, $B$ cannot receive and transmit at the same time, that is, $B$ is operating in a half-duplex mode. If $A$ were to send $y$ to $B$ while $B$ is sending $x$ to $C$, $B$ just ignore the transmission from $A$: $B$'s transmitter is active and so its receiver must be in-operational till $B$ is done talking to $C$.


The Noise Figure of a passive attenuator is the attenuation value of the attenuator. This is because noise figure is by definition the difference in the SNR of the output and the SNR of the input, with the input terminated with the characteristic impedance of the system (typically 50 ohms).  Using this definition, the noise figure of an amplifier similarly is the additional increase in noise beyond the gain of the amplifier when the input of the amplifier is terminated.
The ever-present thermal noise floor is a white noise source with a power spectral density given by kT with k as Boltzmann’s Constant and T as temperature. At room temperature in power units of dBm (dB relative to 1 mW) this ends up being a power spectral density of -174 dBm/Hz. In the case of an attenuator, thermal noise is the only noise source given it is a passive device. If the attenuator is matched to 50 ohms, and with the input tied to a 50 ohm source resistor only and no other gain at the input, the thermal noise at the input will be -174 dBm/Hz (kTB) at room temperature.
The thermal noise at the output of the 50 ohm attenuator will also still be -174 dBm/Hz! The attenuator cannot somehow attenuate the thermal noise that is present in all parts of the system-- it certainly attenuates any noise that is applied to the input, but in terms of power density, the same level of thermal noise is present at the output as a local independent noise source.
Thus in terms of (SNR out) - (SNR in): Consider a noise-free signal at the input with it's noise floor then limited by thermal noise alone; this sets the input SNR. The signal will certainly be attenuated by the attenuation value, but as we explained above the noise can't get pushed down lower than its theoretical limit (the thermal noise floor). Since the noise is not affected as it is the same value at the output and input in this case, the SNR is therefore changed by the attenuator value, which by definition is the Noise Figure.
Cascaded Noise Figure is then considering what occurs when you place an amplifier in front of an attenuator. Now in this case amplified noise is applied to the input, not the thermal noise floor specifically. This amplified noise like any other waveform can of course be attenuated, down to the limit of the thermal noise floor. Consider these two cases and how the overall Noise Figure compares:
The formula for cascaded Noise Figure is as follows below for two components, expandable for any number of components. I'll give the equation first but then demonstrate with two cases where the explanation above becomes very clear.
$$NF_{cascade} = NF_{N-1} + \frac{NF_N-1}{G_{N-1}}$$
This states that the overall noise figure of N cascaded devices is the noise figure of the first component plus the noise figure of the subsequent components reduced by the cumulative gain of to each components input. The gain (G) and noise figure (NF) are given in power ratio (non-dB) quantities.
Consider this with these two components:
Amplifier with 20 dB gain, and 3 dB Noise Figure
10 dB Attenuator
Case 1: Amplifier followed by Attenuator
According to the formula:
$NF_1= 10^{3/10} = 2$
$NF_2 = 10^{10/10} = 10$
$G_1 = 10^{20/10} = 100$
$NF_{cascade} = 2 + \frac{10-1}{100} = 2.09$
The Noise Figure in dB is $10 \log_{10}(2.09) = 3.2 \text{ dB} $
Notice the attenuator had very little effect due to the amplifier in front of it! This is the motivation for Low Noise Amplifiers (LNA) at the front-end of the receiver with very low loss components in front of it (typically a necessary filter).
Case 2: Attenuator followed by Amplifier
$NF_1= 10^{10/10} = 10$
$NF_2 = 10^{3/10} = 2$
$G_1 = 10^{-10/10} = 0.1$
$NF_{cascade} = 10 + \frac{2-1}{0.1} = 20$
The Noise Figure in dB is $10 \log_{10}(20) = 13 \text{ dB} $
Intuition:
Consider a received signal that has a 30 dB SNR, with the noise limited by thermal noise alone. In case 1, we amplify both the signal and the noise by the amplifier with 20 dB gain. This amplifier has a 3 dB noise figure which means the signal goes up by 20 dB, but the noise goes up by 23 dB! (Due to the added electronics noise that is captured in the Noise Figure metric).  So now the SNR at the output of the amplifier is 17 dB, but the thermal noise is amplified and now much higher (23 dB higher!) than the actual thermal noise floor present at the amplifiers output (we have smothered it with an amplified version of the noise). The attenuator can now push both amplified waveforms (signal and noise components down by it's attenuation or 10 dB. Thus the signal that we increased by 20 dB is now reduced by 10 dB, and the noise that we increased by 23 dB is reduced by 10 dB. However the underlying thermal noise is now closer to the amplified noise level, 13 dB away, contributing a small amount to total noise (explaining the additional 0.2 dB in the 3.2 dB Noise Figure result).
Notice the power summation of the amplified thermal noise with the thermal noise 13 dB lower is a total noise increase of 0.2 dB (it all adds up, demystifying cascaded noise figure!):
$$10 \log_{10}(1+10^{-13/10}) = 0.2 \text{ dB}$$


Compute filter output with given values of $n$

Your impulse response is $h[n]=\delta [n] - \delta [n-1]$, so it basically takes a copy of the input and then substracts to it a delayed version of the former. When $x[n]=u[n]$, the output signal is:

$$y[n] = u[n] - u[n-1]$$

This means that we have a step function (i.e. it is equal to $1$ for $n \geq 0$ and $0$ otherwise), and we substract to it a shifted step (i.e. it is equal to $1$ for $n \geq 1$ and $0$ otherwise). Therefore, $\forall n \geq 1$ we are cancelling out the two signals ($1-1=0$). This leaves us with only one not-zero value, which happens to be at the origin. So the output turns out to be

$$y[n]=\left\{
 \begin{array}{ll}
  1  & \mbox{if } n = 0 \\
  0 & \mbox{otherwise }
 \end{array}
\right. \implies y[n]=\delta[n]$$


Why is $x$ not considered a primitive polynomial while being considered an irreducible polynomial?

From the Wikipedia page irreducible polynomials: "a non-constant polynomial $f$ in $F[x]$ is said to be irreducible over [the field] $F$ if it is not the product of two polynomials of positive degree."

From this definition, it should be clear that $x$ cannot be viewed as the product of two polynomials of positive degree. Take one polynomial to be $x$ then there is no other polynomial (with degree greater than zero). Irreducible polynomials act like prime numbers in the field in that other polynomials can be factored into (powers of) the irreducible polynomials.

A primitive polynomial has a different definition. "In field theory, a branch of mathematics, a primitive polynomial is the minimal polynomial of a primitive element of the finite extension field $GF(p^m)$." Notice that irreducible polynomials can be defined over any field, but primitive polynomials only apply to the extension of a finite field. These polynomials are minimal polynomials that generate all of the elements of the extension field.

From this definition, we require that $x$ generates all of the (nonzero) elements of the field. It does in this case. We also require that it is a minimal polynomial which means it has a root which is a primitive element (generates all of the nonzero elements of the extension field). It fails this test as the root of $x$ is $0$ and $0$ is not a primitive element of any field since $0^n = 0$.


Lower rolloff means less excess-bandwidth and therefore longer time span.
That is, you will need a longer filter to accurately reproduce the ideal RRC filter response (increases implementation complexity)
Longer filter means more accurate RRC filter and potentially lower ISI at the receiver side (depends on your target SNR)
Longer filter means less noise out of band transmitted (infinite length filter would have zero power outside symbol-rate*(1+rolloff)

In practice, due to the above constraints, RRC filters are usually not designed by just truncating the ideal RRC impulse response. Either you multiply the impulse response with a window function or you directly optimize the filter coefficients for error and rejection.


A cheap and simple way is to differentiate the signal numerically and threshold the magnitude of differencate to identify the abrupt changes.While deciding the threshold you should take into account the dynamic range of the signal.


A few reasons why this is a terrible idea:

Unless your frequency is an exact multiple of the FFT window size, it'll be spread over adjacent FFT bins.
You're mentioning Morse code, so your carrier is modulated. This is not a constant tone, it has a changing envelope, and thus, its spectrum is convolved by the spectrum of the envelope. As a result, you won't see just a spectral peak at the carrier frequency. Instead, it'll be spread over a frequency band proportional to the symbol rate.
Selectively "editing" FFT coefficients is equivalent to multiplying the spectrum by a function with sharp edges, and thus is equivalent to convolving in the time domain by a window with lots of ripples - causing ringing (More on this here).
Unless you want something ridiculously selective, a bandpass/peaking filter is cheaper in terms of computational cost. Moreover, using the FFT/IFFT approach on long signals (or in realtime) is possible only through the overlap-add method, which is less straightforward to implement than a filter.

Finally, there's no reason why you would want to leave the imaginary part unmodified.


The SDR (or any general digital signal processing system) takes the received RF signal, and downconverts it from the carrier frequency to the baseband.

Now, the real bandpass signal from the antenna does not necessarily have a symmetric spectrum around the carrier frequency, but it can be arbitrary. If the downconverter now shifts the spectrum to the center frequency, the corresponding time-domain signal becomes complex. So, the I and Q samples you get from the SDR are the real and imaginary part of the complex baseband signal, that corresponds to your real passband signal around the carrier frequency.

More details can e.g. be found at the wikipedia site for digital down conversion.

To answer your question:

I/Q representation does not correspond to different sampling points of the signal. Instead, it corresponds to the real and imaginary part of the digital complex-valued baseband signal. These parts are obtained by separately multiplying the RF signal with a sine and a cosine and sampling both streams after low-pass filtering.

Sampling with double frequency can yield the same information as I/Q. It would be neccessary to homodyne the signal to $f_s/4$ to get all information that would have been in the baseband IQ signal to be in the passband signal at $f_s/4$ (where $f_s$ is the sampling rate).


What you did should work, up to a constant, as mentioned by Dan Boschen in a comment. The Hilbert transform will remove any DC component, so the output of the second Hilbert transform should equal the (negative of the) input, up to a constant, and, of course, up to numerical accuracy.

The following example (in Octave) shows this:

t = 0:1024;
x = sin(pi/7*t)+3*cos(pi/5*t);    % some signal
xh = imag( hilbert(x) );
xhh = imag( hilbert(xh) );
err = x - mean(x) + xhh;
max( abs(err) )                   % 6.2172e-15


To stretch the spectrum, you will also need to do some high quality (Sinc kernel) interpolation to fill in any bins that don't have an integer submultiple source.


What is difference between along the edge and across the edge

This is illustrated for a step edge. If you look at gray-scale values, one often considers that a 1D edge is (relatively) irregular across the edge (from gray to black), and regular along the edge (parallel to the border, the luminance value does not vary).

You can perceive it better in the 3D view, where the amplitude along with the gray value is represented in 3D.

In more generic cases, at the appropriate scales, the structure tensor is supposed to capture this behavior.


Fir filter not working according to frequency response

Your normalized stop band and pass band edges are:

$$\begin{align}f_{s1}&=0.5/30=0.017\\
f_{p1}&=1/30=0.033\\
f_{p2}&=3/30=0.1\\
f_{s2}&=4/30=0.133\end{align}$$

And the normalized frequency of your sinusoidal input signal is $0.6593/50\approx 0.013$, which is smaller than $f_{s1}$, so it's in the filter's stop band.

Note that you've designed your filter for a certain sampling frequency, and that it is only the normalized band edges that are fixed; the absolute frequencies of the band edges change according to the sampling rate.


Confused about the period of a discretised time-series

On the continuous setting, the period is the minimum $T\ne 0$ such that $\theta(t+T)=\theta(t)$. Once you decide to have $N$ regularly sampled points inside a period, you introduce a sampling period $T_N$ that is directly proportional to the period: $T_N=T/N$. Hence, you can describe your period as a quantity with a time unit (second) or as a unitless number of samples, equivalently.

Now, if for some reason you normalize the time axis with an affine function, the mean will only shift the whole axis (no change for the period), and the scaling only will dilate or shrink the interval between samples.

So here, the effect of  $s(1)$ will turn $T_N$ into $T_N/s(1)$, and also  $T$ into $T/s(1)$.

In your case, you can directly compute $s(1)=T_N\sqrt{\frac{N(N+1)}{12}}$. In this formula, the scale factor depends on $N$ only, and not on the start and end indices.
For $N=100$, the new period should be now a unitless quantity of
$$T/s(1) = \sqrt{\frac{12N}{N+1}}\,.$$


I understand that when we take longer time windows and take the fourier transform we would suffer in the time resolution, But Why if we take a shorter time window we would suffer in frequency resolution?

Imagine that your signal is two different sinusoids that are close in frequency and they both start at phase 0.  If the time window is short their end phase will be almost identical because their frequencies are close.  Thus, they are difficult to distinguish with a short time window.

If the time window is long the end phase will be quite different- eventually they will be opposite phases (180 degrees out of phase), and thus very easy to distinguish.  That is why longer time windows give better frequency resolution.

  why longer time intervals are needed for high frequency signals? and shorter time intervals are needed for low frequency signals? In another word, Why we want good frequency resolution for high frequencies and good time resoution for low frequencies?

You do not need longer time intervals for high frequency signals.  You'll need to clarify what you're trying to get at on #2 if you want a better response.

  If I have a signal contains frequencies from 0 Hz to 50 Hz, How can I define the range of low and high frequencies?

Again, I'm not sure what you mean.  What you define as a "high frequency" or "low frequency" is completely arbitrary.  What you will need to do is sample the signal at more than 100 Hz to make sure that the Nyquist frequency is above 50 Hz.


Reducing sampling rate by a non-integer factor

I can recommend the resample command in MATLAB. It's a combination of up- and down-sampling. To get from 12kHz to a sampling rate of 9kHz, you upsample by 3 and downsample by 4. resample just needs your signal and the two factors for up- and down-sampling.  Remember, when downsampling, to first low-pass filter out everything between 4.5 kHz and 6 kHz. Very little energy should exist above the resulting Nyquist frequency.


Matlab: is a highpass filter what I want, to obtain my required spectrogram? How do I design/apply this? (pic included)

If you don't have any specification for the filter besides the cut-off frequency, I believe the simplest way is to use the fir1 function. Note that the frequencies are normalised to half the sampling rate. For example:

b = fir1(n,20e3/(250e3/2),'high'); %Wn = Cutoff_req/(Fs/2)
figure; freqz(b,1,[],250e3); %visualize the filter if you wish
FilteredData = filter(b,1,audioData);

If you wish to have more control over the specifications (attenuation, ripple, etc.), I suggest using fdatool.


Is there a relatively easy way to detect likely real-time peaks in discrete-time data?

The most straight-forward solution would be to interpolate the discrete values around the maximum you're interested in, and compute the location of the maximum of the interpolating function. I'd start by using a quadratic or cubic polynomial. This is a simple example of using an interpolating quadratic polynomial.


Gabor filters for digital image processing

In short: If a filter have DC component it means that for a constant region of the image, the filter response will depend on the grey value of that region. Constant region with higher values ==> higher response.

One usually wants higher responses only in case of edges/corners and not in the case of bright but constant regions.

As a picture worth 1000 words, here is an example. I took a simple image consisting of two horisontal lines, each one in different background. I took the Gabor and LogGabor filter. LogGabor is a variant that have 0 DC response.

I filtered the image with a Gabor filter (tuned to respond to horizontal signals) and then, the same image, with LogGabor filter (tuned with exactly the same parameters as Gabor filter).

Original image:

Gabor convoluted:

LogGabor convoluted:

Another example of Gabor filter, with a bandwidth of 0.7 octaves.

Now, let's ignore the edges for a moment and focus on the constant regions.
In the second image, you see different responses for the constant regions! Darker region have a lower response than the lighter one. For the LogGabor case (third image) you see that the response is the same regardless of the original intensity.

Note: The response values were scaled in 0-255 interval, so, neutral grey means 0 response, darker means negative values and lighter positive ones. The bandwidth is set to 2 octaves (except for the 4'th image)

EDIT

According to a paper by Boukerroui et al., for bandwidths around 0.7 the DC response is close to the computational errors ($10^{-6}$). The 4'th figure demonstrate this.

My personal advice to you is to use another bandpass filter. I strongly recommend Gauss Derivatives or Difference of Gaussians (in this order) if it is necessary for you to stay in spatial domain. If you don't have this restriction, I personally prefer LogGabor.

From my experience there is no net advantage or disadvantage in terms of end result between logGabor, Gaussian Derivatives, Difference of Gaussians or more exotic Cauchy filters. BUT this might vary, depending on your practical application.

As for the code, here is a snippet of my Gabor/LogGabor implementation.
The code is in Matlab, the convolutions are done in Fourier space. The complexity is $O(n log n)$ and it depends only on the image dimensions, not on the filter dimensions. The code follows the guidelines in Boukerroui et al.. It is developed by me and released under MIT licence (not shown here) so feel free to use it as you see fit.

Img must be a MxNx1 of doubles, (preferably between 0 and 1 but it should work with 0-255 too); P0 is the tuning period (wavelength), that is, $P_0=1/f$, where f is tuning frequency; Orient is the filter orientation, expressed in radians; FBW frequency bandwidth in octaves; ABW the angle bandwidth.

function [ ReConv,ImConv ] = GaborConv( Img, P0,Orient, FBW, ABW,varargin )
 ReConv=Img; ImConv=Img;
 ImgFFt=fft2(Img);
 lg2pi=sqrt(log(2)/pi);
 A=P0*lg2pi*(2^FBW+1)/(2^FBW-1); B=P0*lg2pi*(1/tan(0.5*ABW));
 K=1/(A*B); KF=2*pi/P0;
 [N,M]=size(Img); u0=1/P0*cos(Orient); v0=1/P0*sin(Orient);
 FFilt=zeros([M,N]);
 stx=1/(M-1); sty=1/(N-1);
 [U V]=meshgrid(-0.5:stx:0.5,-0.5:sty:0.5);
 Ur=(U-u0)*cos(Orient)+(V-v0)*sin(Orient);
 Vr=-(U-u0)*sin(Orient)+(V-v0)*cos(Orient);
 FFilt=exp(-pi*(Ur.*Ur.*A.*A+Vr.*Vr.*B.*B)) ;
 FFilt=ifftshift(FFilt);
 normFact=sum(sum((abs(ifft2(FFilt)))));
 RezFreq=FFilt.*ImgFFt/normFact; RezSpace=ifft2(RezFreq);
 ReConv=real(RezSpace); ImConv=imag(RezSpace);
end

Log Gabor method:

function [  ReConv,ImConv ] = LogGaborConv( Img, P0,Orient, FBW, ABW,useCos ).
 ReConv=Img;
 ImConv=Img;
 ImgFFt=fft2(Img);
 %Compute the k_beta according to the desired frequency bandwidth
 kBeta=exp(-1/4*sqrt(2*log(2))*FBW);
 gaussSD=sqrt(ABW^2/(2*log(2))); %compute Gauss SD to ensure that Gauss(ABW,0,SD)==0.5
 lnkB2=2*log(kBeta)^2;
 nc=exp(-1/8*lnkB2/2)*sqrt(-2*sqrt(pi)/(1/P0*log(kBeta)));
 [N M]=size(Img);
 stx=1/(M-1);
 sty=1/(N-1);
 [Ug Vg]=meshgrid(-0.5:stx:0.5,-0.5:sty:0.5);
 Radius=sqrt(Ug.^2+Vg.^2);
 Radius = ifftshift(Radius);
 Radius(1,1)=1;
 Radius = fftshift(Radius);
 %prepare a low pass filter, as proposed by Kovesi
 lowpass=1./(1+(Radius./0.45).^(30));
 %Compute the bandpass filter
 logRad=nc*exp(-log(Radius*P0).^2/lnkB2);
 RadialBP=logRad.*lowpass;
 %Compute the angle component
 Theta=atan2(Vg,Ug);
 sintheta = sin(Theta);
 costheta = cos(Theta);
 ds = sintheta * cos(Orient) - costheta * sin(Orient);
 dc = costheta * cos(Orient) + sintheta * sin(Orient);
 dtheta = abs(atan2(ds,dc));
 if ~exist('useCos','var')
     useCos=0;
 end
 if useCos==0
     %Use a gauss angle
     AngFilter=exp(-dtheta.^2/(2*gaussSD^2));
 else
     %Use a cosine angle
     AngFilter=cos(dtheta).^2.*heaviside(dtheta./pi+0.5).*heaviside(0.5-dtheta./pi);
 end
 logGaborFilt=RadialBP.*AngFilter;
 logGaborFilt=ifftshift(logGaborFilt);
 logGaborFilt(1,1)=0;
 normFact=sum(sum((abs(ifft2(logGaborFilt)))));
 RezFreq=logGaborFilt.*ImgFFt/normFact;
 RezSpace=ifft2(RezFreq);
 ReConv=real(RezSpace);
 ImConv=imag(RezSpace);
end

As for how the DC response varies with various parameters and why, a comprehensive discussion is in Boukerroui et al., equations (19) and (20) at pages 75, 76. Take a peek at Appendix C (pag 101). I would like to redirect you to the paper rather than just copy/paste the relevant paragraphs.

Hope it helps!


How to implement a digital sawtooth oscillator?

Refer to think link, and use the truncated phase word directly (no look-up-table needed):

Numerically Controlled Oscillator (NCO) for phasor implementation?


We can assume
$g=i(x,y)*f(x,y)+n(x,y)$, where $g$ is the above images, $i$ is the original image , $f$ is a degradation function and $n$ is image noise. If you find $f$ and deconvolve it with $g$, you'll get $i$ plus some noise. Note, simple deconvolution would not work, because the noise ($n$) in the process is amplified. The Wiener filter can used to restore images that have been degraded by a known degradation function.Here I am not sure what is degradation function.
I suggest you read the following material as well:
Lecture notes of Prof. William Hoff


Impulse response of forward difference cascaded with one sample delay

For any sequence $x[n]$ you have
$$x[n]\star\delta[n-n_0]=x[n-n_0]\tag{1}$$
because
$$x[n]\star\delta[n-n_0]=\sum_{k=-\infty}^{\infty}x[k]\delta[n-n_0-k]\tag{2}$$
Note that $\delta[n-n_0-k]$ is only non-zero for $k=n-n_0$, hence the result $(1)$.
Now use $x[n]=\delta[n+1]-\delta[n]$ and $n_0=1$.


Disparity in phase calculation between STFT and Hilbert transform

Without seeing your data, it's hard to tell whether you have contributions of closer frequencies or not.  That is certainly possible.  So, my suggestion to you would be this: try out your two different methods on pure sine waves where you know that there are no closer frequencies.  In other words, generate some dummy data which is just one single frequency.  If your two methods give the exact same answer on that data, then maybe it is in fact the contributions of closer frequencies that are messing up your original calculation.  But if the two methods do not give the exact same answer, then you likely have an error in your method.


It's a version of Parseval's Theorem, which can be easily proved by noting that

$$\int_{-\infty}^{\infty}s(t)r^*(t)dt=\mathcal{F}\{s(t)r^*(t)\}\big|_{f=0}\tag{1}$$

where $\mathcal{F}$ denotes the Fourier transform. Note that in general the integration limits should be $-\infty$ and $\infty$ for both integrals. I assume your notation refers to band-limited functions.

From $(1)$, the proof is quite short:

$$\mathcal{F}\{s(t)r^*(t)\}=S(f)\star R^*(-f)=\int_{-\infty}^{\infty}S(\nu)R^*(\nu -f)d\nu\tag{2}$$

where $\star$ denotes convolution. So

$$\mathcal{F}\{s(t)r^*(t)\}\big|_{f=0}=\int_{-\infty}^{\infty}S(\nu)R^*(\nu)d\nu\tag{3}$$

which establishes the equation in your question.


How to demodulation BPSK in real case(through a channel)?

There will inevitably be frequency offsets and time offsets in the OP’s receiver, given the receiver is working on a completely independent clock, so cannot have the exact same carrier frequency (and in the case of mobile stations there is also Doppler offset). The typical BPSK receiver will use carrier and timing recovery loops to extract the estimates of the carrier and time offset from the signal itself.
A very simple carrier recovery for BPSK is done by simply squaring the received signal, which will produce a dominant synchronized tone at twice the carrier frequency. On its own this would usually be too noisy to use as the carrier, so is typically locked to a clean local oscillator using a PLL and then dividing the resulting frequency by two, which is then multiplied by the received signal to demodulate the BPSK signal.  If one has the flexibility of using a carrier that is an integer multiple of tour  symbol clock, then timing recovery can be done simply with additional frequency dividers; otherwise the Gardner Timing Loop described in the links below can be considered (an error discriminator is needed such as the Gardner Timing Error Detector; this error then gets accumulated and the accumulated error weighted by the desired loop gain controls the sampling location for the timing clock resulting in true recovered data).
The recovered carrier will have 0/180° phase ambiguity so further information to resolve this is encoded in the data itself. Data can be encoded with a known preamble that has many data transitions to allow for resolving this ambiguity and also to answer the OP’a multipath question, help train a channel equalizer. An example preamble could be a 11 bit barker code (10110111000) followed by (1010101010), or depending on acquisition SNR required could be a longer series of these or similar codes. The barker code has the excellent autocorrelation property that it is not highly correlated with shifted versions of itself so can help establish the correct start of a data packet. This autocorrelation property (ideally an impulse at $\tau$ = 0 ) is also ideal for training an equalizer since it is approximately white in frequency: An equalizer can only compensate for the frequencies it has measurements for. The 101010...sequence is very friendly to timing recovery given the transitions on every bit.
In addition the signal level will likely be variable so one would likely want to have an AGC (automatic gain control) loop as well.
Other posts where I provide more details on carrier and timing recovery and equalization all applicable to BPSK are summarized below:
Carrier Recovery:
Phase synchronization in BPSK
FFT-based coarse carrier recovery for QPSK
High modulation index PSK - carrier recovery
Recovering signal for psk
Timing Recovery:
High modulation index PSK - carrier recovery
Location of Matched Filter
Isn't Gardner's algorithm and Early-Late gate the same thing?
Channel Equalization:
Compensating Loudspeaker frequency response in an audio signal
non linear equalizer vs linear equalizer
Does zero-forcing equalizer need known channel impulse response?
The benefits of a fractionally spaced equalizer
Recursive Least Square Adaptive Linear Equalizer


Avoiding identifying objects within previously identified objects in MATLAB

There are two alternatives to this.

Adjust your opening parameter so that you end up with a fully connected shape. Say, instead of 15 you can use 25.
You can find the areas of all intersecting polygons using e.g. this. If the area of an intersection is higher than a given threshold, discard the polygon.

After finding overlapping polygons, you could consider averaging the hypothesis so that you end up with a more correct shape pose.


Solve equations with conjugates multiplications

can we solve x2 and x1 if y2 and y2 are known?

Yes, provided that you mean $y_1$ and $y_2$ are known. If you have two linear equations your can (in most cases) solve for two unknown variables.

If not, can we get the value of y2 based on the equation (1),

No. $y_2$ doesn't show up in equation (1) at all, so you can't get it from there.

which is very similar to equation (1)

Similarity has very little to do with it. You can always solve for two unknown variables unless the two equations are linear dependent. The two that you have are not.


How to convert A/D data to the needed representation

Two methods:

(1) theoretical - get the specs for your load cell, analogue hardware, A-D, etc, and do the math

(2) practical - calibrate your system by applying a series of known accurate loads, logging the A-D output and then derive a calibration offset and gain (assuming everything is supposed to be linear)

Ideally you should do both, i.e. do (1) first and then do (2) to calibrate out any offset or gain errors in your system.