This depends on your modulation choice and approach to carrier recovery. But consider this for a specific case of a raised cosine pulse shaped signal:

Observe (using the eye diagram would be simple) the significant symbol to symbol jitter everywhere except in close proximity to the symbol decisions — for this case I would opt to use 1 sample per symbol (the decisions) to estimate symbol to symbol phase since those samples have the least jitter (this would be the samples after all equalization, timing recovery and matched filtering is complete, just prior to decision where the complex phase can still be measured between two successive symbols).

This approach works well if you are using the Gardner TED for symbol clock recovery since it can function well with relatively large carrier offsets.

You can evaluate your specific waveform with the same considerations to determine what approach would have the minimum estimation error (and not the error on any given sample but the rms error within your carrier tracking loop bandwidth).

Here is an example shown below for a raised-cosine filtered waveform. Typically the filtering is split between the transmitter and receiver with Root-Raised Cosine (RRC) filters such that sufficient pulse shaping for reducing out of band emissions is provided while allowing for a matched filter for optimum SNR (against white noise conditions) in the receiver. In the plots below you see the signal as received (if noise-free) and then after it is passed through the RRC matched filter.  Notice the wider variation in the zero crossings after the RRC filter in the plot on the right compared to the plot on the left. It is the zero crossing locations that are used, on average, to determine the data clock. This is why, without additional processing to minimize the zero crossing variation (which can be done), it is better to use the the waveform on the left for timing recover. However notice the sharp points at the correct sampling locations for symbol decision on the plot on the right compared to the one on the right-- these are the best samples (at one sample per symbol using the best derived sample from timing recovery) to measure phase rotation for carrier recovery given the Gardner TED can operate with significant carrier offset.


Proving that a product of matrices invertible

First, this question is probably better for MATH.SE, but I'll give it a shot.

It's been a long long time since I did this stuff.

If $N > M$:

1) $C^*$ has rank M.

2) $R_x^{-1}$ has rank M. (or it wouldn't exist as an inverse)

3) Therefore  $C^*R_x^{-1}$ has rank M, since $R_x^{-1}$ is a full rank square matrix.

4) $C$ has rank M

5) $C^*R_x^{-1}C$ has rank less than or equal to M and is NxN.

6) $C^*R_x^{-1}C$ is a not a full rank square matrix and is thus not invertible.

If $N \le M$:

1) $C^*$ has rank N.

2) $R_x^{-1}$ has rank M. (or it wouldn't exist as an inverse)

3) Therefore  $C^*R_x^{-1}$ has rank N.

4) $C$ has rank N

5) $C^*R_x^{-1}C$ has rank less than or equal to N and is NxN.

6) $C^*R_x^{-1}C$ could be full rank square matrix and is thus could be invertible.

At this point, I am not sure what is required for $C^*R_x^{-1}C$ to have rank N.  I have narrowed the results, but not fully answered your question.

Gosh, I hope I have this right.

https://math.stackexchange.com/questions/1524444/connection-between-rank-and-positive-definiteness

https://math.stackexchange.com/questions/272049/rank-of-matrix-ab-when-a-and-b-have-full-rank

Starting with another disclaimer:  Ordinarily I don't answer questions here unless I am fairly rock solid in my understanding of the subject matter.  In this case, as the edit history shows, I am muddling about quite a bit.

I had to look up the Frobenius inequality and I understand how you get $ M \leq \text{rank}(C^* R_x^{-1} C) $ from it in the $N \ge M$ case.  I still don't see how it could be greater than $M$.

Oriol B inserted:

  Regarding your question about the Frobenius inequality, I derived it as follows:
  $\text{rank}(C^* R_x^{-1} C) \geq \text{rank}(C^*R_x^{-1}) + \text{rank}(R_x^{-1}C) - \text{rank}(R_x^{-1})=\text{rank}(C^*R_x^{-1}) + \text{rank}(R_x^{-1}C) - M $. And by now using Sylvester's rank inequality, we get $\text{rank}(C^*R_x^{-1}C) \geq \text{rank}(C^*R_x^{-1}) + \text{rank}(R_x^{-1}C) - M \geq rank(C^*) + \text{rank}(R_x^{-1}) - M + \text{rank}(R_x^{-1}) + \text{rank}(C) - M - M = \text{rank}(C^*) + \text{rank}(C) - M = M $
  and the last inequality follows by the assumption.

I do feel confident in my $N>M$ argument, but thanks to this addition, you can change #5 to saying the rank of $C^* R_x^{-1} C$ is $M$.  However, if $ N > M $, then the rank is still smaller than the dimension and it is still not invertible.  It should be with $N=M$.

I'm still mucking with the $N \le M$ case.  Using Frobenius you get $ \text{rank}(C^* R_x^{-1} C) \ge 2N - M $.  For $N=M$ it looks like the result will have rank N and be NxN so it will be invertible.

Since the $N<M$ case is of no interest to you I am stopping for a while on this.

https://artofproblemsolving.com/community/c364309h1480887_rank_inequalities_and_some_consequences

Another update:

I think I found the first flaw in your argument.  I mean, there has to be one in order for mine to be correct.  ;-)

$C$ is $M$x$N$ and $y$ is $N$x$1$ so $Cy$ is $M$x$1$.

$Cy$ can be interpreted as the linear combination of $N$ $M$ length vectors (the columns of C).  If $N>M$, your assertion that Cy can't be zero is false.

I'm not sure where your misconception arises.  Here is a concrete counter example:

Suppose $M=2,N=3$

$$
\begin{bmatrix}
1 & 0 & 1 \\
0 & 1 & 1
\end{bmatrix}
\begin{bmatrix}
1 \\
1 \\
-1
\end{bmatrix}
=
\begin{bmatrix}
0 \\
0
\end{bmatrix}
$$

$$ Cy = 0 $$

Clearly, $y \ne 0$ yet $Cy$ is.  The rank of $C$ is 2.

Response to OriolB's "Of course your assumption ..." comment:

By "your assumption" you mean my second case of $N \le M$.  I didn't get that, sorry.

It seems that you have removed the "less than or" part of step 5, thus proving it invertible.  I had left it hanging since you had said you were only interested in $N>M$ which I proved was not invertible in the first case.

For the $N>M$ case, perhaps a pseudo-inverse will suffice for your greater purposes.


20Vrms is the maximum voltage you can apply across the transmitting transducer without the risk of immediately damaging it.
The amplitude of sound it produces is determined by the driving voltage. The transmitter is characterized at 10Vrms, so about 28Vp-p assuming a sine wave, probably where you would prefer to use it for reliability and long life.


What does imaginary Component of a bin's phase in the Fourier Transform mean?

It simple trigonometry.  A vector (FFT result bin) can have a magnitude and a phase.  Or you can specify the same vector identically with X,Y coordinates.  The complex representation simply uses X for the real part and Y for the imaginary component.  Without an imaginary component, there’s no place to put the Y, and your vector would be under specified.

The phase ends up being atan2(I,Q) == atan2(Y,X) == atan2(im(z),re(z))

As for why use a complex representation rather than 2 strictly real quantities?  The complex math is more elegant, and requires less writing on the chalkboard to demonstrate similar things.


Find frequency/phase response of a microphone

I would say skip MATLAB and learn more about the published frequency responses by looking up your mics here: http://recordinghacks.com/microphones. There are hundreds of mic models with responses as such:


This would appear to be a bug in the freqz implementation. I didn't dig enough to figure out why but even this doesn't work:

freqz(SOSb(1:2,:));

Even though angle(freqz(SOSb(1:2,:))) returns a legitimate result.

But I found this as a workaround:

[h,f,s]=freqz(SOSb);
freqzplot(h,f,s);

Unfortunately freqzplot is being obsoleted so you might have to resort to:

plot(angle(freqz(SOSb)));


How to change speed of audio samples without changing pitch?

This question (about "time-scaling" audio) is closely related to pitch shifting, which is time-scaling combined with resampling.  But changing the speed without changing pitch is only time-scaling, so there is no resampling involved (contrary to what thomas has suggested).

There are frequency-domain methods (phase-vocoder and sinusoidal modeling) that can change speed for an orchestral mix (or some other broadbanded non-periodic sound) without glitches, but they can become a little "phasey" in their stretched sound.

If it's monophonic, like a solo or a single note or string or voice, a less expensive time-domain approach is sufficient and can sound very good, under this monophonic condition.  If it's a single note or tone, the signal is quasiperiodic.  That means, in the vicinity of some time, $t_0$, that

$$ x(t) \approx x(t - \tau(t_0)) \qquad \text{for } t \approx t_0 $$

So $\tau(t_0) > 0$ is the estimated period of the quasiperiodic waveform in the neighborhood of $t_0$.  This estimated period $\tau(t_0)$ is estimated by finding the value of $\tau$ that minimizes something like

$$ Q_x(\tau, t_0) = \int \Big|x(t) - x(t-\tau)\Big|^2 w(t-t_0) \, dt $$

where the window $w(t-t_0) \ge 0$ is centered around the time $t_0$.  If there are multiple candidates of $\tau$ that minimize $Q_x(\tau, t_0)$, usually it is best to pick the smallest $\tau$ in which $Q_x(\tau, t_0)$ is small.

Now imagine copying your note $x(t)$ and offsetting (delaying) the copy in time by the amount $\tau(t_0)$.  That delayed copy would be $x(t - \tau(t_0))$.  Then you have two identical waveforms, except for the offset of $\tau(t_0)$ but around the time $t_0$ the two waveforms will look almost the same because the offset is exactly one period, based on the estimate of the period around time $t_0$.

Then you can cross-fade (using a Smoothstep function $S(t)$) from the original to the offset copy, the cross-fade will not have any nasty cancellation or "destructive interference" because the waveforms are lined up.  And the result will be the same note, but one period longer.  (Longer by $\tau(t_0)$ seconds.)

$$ y(t) = \Big(1-S\big(\tfrac{t-t_0}{\tau(t_o)}\big) \Big) \, x(t) \ + \ S\big(\tfrac{t-t_0}{\tau(t_o)}\big) \, x(t - \tau(t_0)) $$

(That won't stretch it long enough to notice a difference, but if you repeat this operation many times per second, you can stretch a patch of sound by even a factor of two.)


Most modern speech recognition systems are based on Hidden Markov Models. These are probabilistic models that require a previous training.
So here is simple comparison:

By using phonemes you have a reduced number of possibilities (around 50 depending on the language), you get a lot of data to train each phoneme model, by combining several models you can have word models.
By using word models since the beginning, you would need to have thousands of models to train and probably not enough data.

If you are working on a very simple recognition system, with a very limited set of words, for example, recognition of numbers, then maybe you can use words models that don't require HMMs.


Compute phase coherence between two eeg signals

Perhaps the built-in 'mscohere' and 'cpsd' functions may help you.

The mscohere function returns a value between 0 and 1 that measures the correlation between the signals, and the phase delay can be computed using the cpsd function, as per this example from the Mathworks website. It suggests that the relative phase between the correlated components can be estimated with the cross-spectrum phase.

https://uk.mathworks.com/help/signal/examples/measuring-signal-similarities.html#d119e1563

Taken from the Mathworks page above

[Cxy,f] = mscohere(sig1,sig2,[],[],[],Fs);
Pxy     = cpsd(sig1,sig2,[],[],[],Fs);
phase   = -angle(Pxy)/pi*180;


A matched filter is matched to the received pulse shape. If the matched filter is implemented in the baseband, i.e., after down-conversion, the received pulse shape can be complex-valued, e.g., if the channel frequency response is not symmetric with respect to the carrier frequency. In that case the matched filter is a complex-valued low-pass filter matched to the received (complex low-pass) pulse shape. The frequency response of such a filter is centered around DC but it is not symmetric with respect to DC.


Discrete Fourier transform doesn't actually decompose a signal into sum of sinusoids it was composed from. DFT projects signal onto a specific set of discrete sinusoids known as Fourier basis. Sinusoids used to compose the signal may or may not coincide with the DFT basis. If they are not, they are being decomposed into a number of DFT basis elements.

There's another way of how to look at it. DFT deals with periodic functions. If you plot 3-4 periods of your signal interval, you'll see it doesn't look like a sum of pure sinusoids and has discontinuities - if sinusoids it was composed from don't coincide with a Fourier basis.


How to convert time-domain signal to complex envelope?

These are the main ideas:

Consider a receiver that picks up a signal $r(t)$. This signal has bandwidth $W$ and is centered on carrier frequency $f_c$.
Using the Hilbert transform, eliminate the negative frequencies of $r(t)$. The resulting signal, $r_+(t)$, is called an analytic signal.
Now, downconvert the analytic signal using a complex exponential, so that the spectrum is centered around 0. This signal, called $\tilde r(t)$, is the complex envelope of $r(t)$. Its spectrum goes from $-W/2$ to $W/2$, and its bandwidth is $B=W/2$.
The complex envelope is almost always complex, but if the spectrum of $r(t)$ is symmetrical around $f_c$, the CE is real. Examples of this are AM DSB and BPSK.

In the transmitter you would do the opposite:

Start with a complex envelope $\tilde s(t)$ that represents the information you want to transmit. Its spectrum should go from $-W/2$ to $W/2$. An example would be a QAM signal with complex symbols.
Upconvert the CE to frequency $f_c$, multiplying by a complex exponential. The result is an analytic signal with only positive frequencies. The signal will be centered around $f_c$ and cover the frequency range from $f_c-W/2$ to $f_c+W/2$.
Take the real part of the analytic signal to convert it to a real signal that can be physically generated and transmitted.

All of these steps can be accomplished in Matlab. Be sure to read the documentation for the hilbert command. I can't comment on the corresponding Simulink blocks, since I don't use them.

Also note that it is very likely that you don't have to do all this. If you can perform all your simulations using the complex envelope, you'll save a lot of computing time. The reason is that the transmitted and received signals need to be sampled at least at a rate $2(f_c+W/2)$, whereas the complex envelope can be sampled at $W$ real samples per second -- a huge difference for large $f_c$.


Why are some JPEG files of the same scene much larger than others?

My guess is on autofocus or aperture briefly changing in a manner where the resulting image contains more high-frequency elements.

For example, if the focus moves from a smooth object to a textured object (like the smooth hippo to the cloth drape) in a manner where the latter's detail texture surfaces, JPEG will tend to take quite a larger size.

As someone else already stated: it seems like a good idea to check the image EXIF data for changes in core parameters like aperture and focus distance.  For such a striking difference in image size, it is quite likely that some fundamental parameter is different in the camera's opinion.


What Kind of Features Can I Extract from a Signal

Some Features:

Mean.
Variance.
Skewness.
Kurtosis.
Dominant 3 frequencies in the DFT.
Energy of the 3 dominant frequencies.
Max Value.
Min Value.
Median.
Total Variation.

Usually I'd compute them in running windows.
Another great information is the Histogram of the Derivative.
Or just all the above of the Derivative.


Get an input signal from LTI system output, find impulse response

Ok.  You are touching on the topic of deconvolution.  You can do this most easily by transforming the data into the frequency domain, doing some basic math and then converting back to the time domain.

For example for case 1, if you transform the signals you have
$ H(f) = Y(f)/X(f) $ so $h[n] = IDFT[Y(f)/X(f)]$

You can do similar manipulations to handle case 2.

I'm not sure case 3 is possible.  I think you need more information to produce x[n]

Here is a site that discusses deconvolution:  http://terpconnect.umd.edu/~toh/spectrum/Deconvolution.html

This is a VERY simplified treatment of this question.  I'm not sure what level of detail is being requested in the posted question.


DFT for audio classification - signals of different lengths

As you explicitly mention, that you are aware of a loss of information and are willing to accept that, a possible solution will be modulo-N reduction or wrapping, as described by Orfanidis in the freely available book "Introduction to Signal Processing" on pages 489ff.

In modulo-$N$ reduction, a vector $x$ with length $L > N$, where $N$ is the DFT length, is divided into non-overlapping blocks of length $N$. These blocks are summed up, resulting in one vector $\tilde{x}$ of length $N$ .

Figure taken from Orfanidis (linked above), Fig. 9.5.1 on page 489.

You then apply the DFT on $\tilde{x}$, which results in $\tilde{X} = \text{DFT}\{ \tilde{x} \}$ with length $N$. Orfanidis proves that $\tilde{X}$ is exactly equal to the $N$-point DFT of $x$, so is equal to $X$.

In your case, that would mean you could find a minimal length $N$, which is guaranteed for all input signals. You then apply a modulo-$N$ reduction on all your input signals and compare these length-$N$ DFTs. Of course, a higher resolution would be possible for longer signals, and you just ignore this by using the modulo-$N$ reduction method.

Similarly, you could just zero-pad all signals to a maximal length $M$. All your signals will then have a high frequency resolution defined by $M$. Important point: you do have that computational frequency resolution, but the physical frequency resolution is still determined by the length of your signals $L$. More details on this are given in any decent signal processing class or on pages 482ff in Orfanidis.


Is cutting EQ scientifically different than boosting EQ?

There is a very compelling reason for preferring cutting over boosting. Assuming you are using the typical bell-response parametric EQ bands, then you will see a significant phase distortion around the resonant frequency, i.e. the frequency you're boosting or attenuating. (There is an exception to this rule: Linear phase EQs that are often used for mastering do not show this phenomenon. Those only exist as digital implementations however and have other problems.)

The most noticable effect of the phase distortion is that it comes with a certain group delay around the resonant frequency, and especially narrow changes can result in audible post-ringing and temporal distortion of transients.

Now if you're boosting what you want to hear more, you introduce this kind of problem exactly in the band that you want to give more prominence. Cutting also introduces the same issues, but you usually cut what you don't want to hear as prominently. So cutting will move the phase artefacts to the sonically less important frequency range while boosting will place it right where you don't want it.

That's the theory at least. In practice, the rule can be relaxed in most cases. The negative effects of resonant phase distortion greatly depend on the derivative of the phase response with respect to the frequency. That means that resonances with very small bandwidth will result in greater issues than wide band resonances. And this is absolute bandwidth, not relative bandwidth. This means for a typical EQ band that has a bandwidth proportional to the center frequency ("constant-Q filter"), the bandwidth automatically gets larger for higher frequencies.

So a bell EQ at a few kHz will have a large bandwidth, even if it has a rather high resonance. Therefore cutting vs boosting doesn't make a huge difference for higher frequencies. It does however get relevant for frequencies in the lower mid and the bass range, especially when you're using high Q-factors.


Regression vector size for prediction, reconstruction and filtration with adaptive filters

Suppose that you are adapting $w$ to minimize $\text{E}(y[n]-w[n]*u[n])^2$ where
$$y[n]=h[n]*u[n]+\nu[n]$$
$y[n]$ and $u[n]$ are known and $\nu[n]$ is an additive noise component.
With a long enough FIR filter you can model any linear system $h$ with any accuracy. You can opt for a shorter filter length in expense of limited modelling accuracy. A key parameter is the sampling rate of the signals. If the signals are over sampled unnecessarily then the required length for $w$ increases.

Suppose that we settle for a certain modelling accuracy and we want to find $w$. If all the samples of $u[n]$ are available and $h$ remains constant you can find an estimate of $h$ namely $w$ using a least square estimation. This estimation will have a limited accuracy. But it is the best you can do. This error is referred to as the Least square error which here I (a bit wrongly) consider the same as Mean Square Error (MSE error).

Since $h$ may change, we need to continuously estimate $h$, and that is why we are using the adaptive filters. Let's say that we use the LMS adaptation to find $w$.
Now on top of the MSE error we have extra error owing to the way that the LMS adaptation is performed. This error is called  misadjustmentis and it is caused by the limited window where $u[n]$ is used in the adaptation. If the LMS step size $\mu$ approaches zero we are including all sample of $u[n]$ in estimation, but then there is no adaptation. It can be proved (Take a look at any adaptive filters book) that misadjustment is proportional to

The filter length
The step size

The other issue is the tracking behavior of the filter. Since $h$ is changing we like to adapt the filter as fast as possible. To increase the adaption speed we can increase $\mu$. To keep the filter stable we should have
$\mu<\frac{1}{\lambda_{\text{max}}}$
where $\lambda_i$ are the eigenvalues of the covariance matrix of $u[n]$.
The situation is worse if the sample of $u[n]$ are correlated (apparantly your scenario). In this case the covariance matrix has different eigen values and the adaptive filter converges with different modes where the time constant $\tau_i$ for each mode will be
$$\tau_i=\frac{1}{2 \mu \lambda_i}>\frac{\lambda_{\text{max}}}{2  \lambda_i}$$
The ratio of largest to smallest eigenvalue $\frac{\lambda_{\text{max}}}{\lambda_{\text{min}}}$ increases as the samples of $u[n]$ become more correlated and the LMS filter becomes slower. Unfortunately the speed cannot be increased idefinitely by increasing $\mu$. Since $\mu$ is limited from above by $\frac{1}{\lambda_{\text{max}}}$.

Back to your case:
It seems that your data is highly correlated that is why removing 14/15 of samples does not change the performance that much.
Throwing away samples of $u[n]$ can help in your case for the following reasons:

You are decreasing the length of the part that you are adapting. This can decrease the misadjustment.
You are whitening your data by down sampling it. This decreases the eigen ratio and improves the tracking behavior of the filter.
Another issue with highly correlated $u[n]$ is the potential divergence of filter at frequencies where $u[n]$ has almost no component. If $u(e^{j\omega_i})=0$ then $h(e^{j\omega_i})$ can go to infinity if the filter length is infinite. The shorter the filter the least likely that this occurs. This problem however can be solved by introducing leakage.

The negative impact of throwing way samples of $u[n]$ is that your are imposing  a certain assumption on $w$. By using only $[u(k), u(k-15), u(k-30), u(k-45), u(k-60)]$  you are assuming that $w$ takes a form like:
$$w=[w_0~...~0~w_{15}~0~...0~w_{30}~...~w_{60}]$$
This limits the ability of $w$ to model an arbitrary $h$.


Oscillations in a first-order discrete-time linear systems

For continuous-time systems, a pole at location $s_0=\sigma_0+j\omega_0$ will create a time-domain contribution of the form

$$e^{s_0t}=e^{\sigma_0t}e^{j\omega_0t}\tag{1}$$

which is a damped oscillation if the pole is in the left half-plane (i.e., $\sigma_0<0$), and if the pole is not on the real axis (i.e., $\omega_0\neq 0)$. For $\omega_0=0$ there is no oscillation. This is the case for a real-valued first-order system, which can only have a real-valued pole.

For discrete-time systems, a pole at $z_0=re^{j\omega_0}$ ($r\ge 0$) causes a sequence of the form

$$z_0^n=r^ne^{j\omega_0n}\tag{2}$$

If the pole is inside the unit circle (i.e., $r<1$) that sequence will decay in amplitude. Note that the sequence will not oscillate only if $\omega_0=0$, i.e., if the pole is on the positive real axis. If the pole is real, but if it is negative, we have $\omega_0=\pi+2\pi k$, i.e., we have an oscillation with maximum frequency, corresponding to an alternating sequence.

In sum, for continuous-time systems if the pole is on the real axis, it will not cause oscillations. For discrete-time systems, there are no oscillations only if the pole is on the positive real axis. A pole on the negative real axis causes oscillations with maximum frequency. Consequently, a discrete-time real-valued first-order system can exhibit oscillations if the pole happens to be on the negative real axis.

For additional insight consider the exact mapping between the $s$-plane and the $z$-plane:

$$z_0=e^{s_0T}\tag{3}$$

where $T$ is the sampling interval. Inversion of $(3)$ gives

$$s_0=\frac{\ln(z_0)}{T}\tag{4}$$

A $z$-plane pole $z_0=-r$ ($r>0$) on the negative real axis maps to a complex valued pole in the $s$-plane:

$$s_0=\frac{\ln(-r)}{T}=\frac{\ln(r)}{T}\pm j\frac{\pi}{T}\tag{5}$$

The frequency $\omega_0=\pi/T$ corresponds to half the sampling frequency, i.e., the maximum frequency of the corresponding discrete-time system.

As a final note, oscillations have nothing to do with whether the system is stable or not (stability in the bounded-input bounded-output sense). For stability, it only matters if the contributions of the poles to the output signal decay or not.


How do you convert angular frequency to $\textrm{Hz}$ from MATLAB's $\tt freqz$ function?

After a little hard thought, I believe I stumbled across the answer. In order to know how to convert, you really need to pay attention to units.

When you use the 'butter' function, it is necessary that you normalize the frequency. By normalizing the frequency, you are essentially removing the time dependency and are only focused on relative values. Commonly in DSP, this frequency is achieved by dividing the analog frequencies by the sampling rate. In this, the units of the normalized frequency are cycles/sample. This was achieved by dividing the analog frequency in hz (cycles/second) by the sample rate (samples/second). Dividing these yields the units of cycles/sample. Notice now that there is no units of seconds, we have thus removed the time dependency.

However, MATLAB is different in the sense that it likes to normalize the analog frequencies in terms of the Nyquist rate as opposed to the sample rate. The Nyquist rate is 1/2 the sample rate. The reason for this is due sampling theorem which states that the sampling rate must be twice that of the true frequency in order to capture the true frequency. For analog frequencies greater than 1/2 the sampling rate, they will be aliased onto frequencies less than 1/2 the sampling rate. Again, just refer to the sampling theorem.

So now, using MATLAB's convention, the real units of the normalized frequency is (1/2 cycles)/sample. Fine and dandy.

After using the 'butter' function to obtain the filter coefficients, we feed them into the 'freqz' function which will create an h and w vector. h and w are the frequency response and angular frequencies respectively. Assuming you know about frequency responses/bode plots, the h vector should make sense.

However, the real question is regarding the w vector. MATLAB states that the w vector is the angular frequency vector. But what are the units of this angular frequency? Well typically when you think of angular frequency in terms of radians, you think of radians/second. However, since we're working with normalized frequencies with removed time dependencies, the units of normalized angular frequencies is radians/sample. If you look at the range of the w vector, you will notice that is goes from 0 radians/ample to pi radians/sample. This make sense due to the fact that pi radians/sample is equivalent to (1/2 cycles)/sample which covers the entire normalized frequency range up until the normalized Nyquist frequency.

So we're now at the point where we have our filter response magnitude 'h', but it's been evaluated at the weird units of radians/sample. We want to get back to our analog frequency equivalent which has units of Hz (Hz = cycles/second).

It's now simply a game of unit analysis of converting radians/sample -> cycles/second.

That's achieved by doing the following operation to the w vector:

w*Fs/(2*pi) -> units of Hz.

TL;DR
freqz computes the digital filter response in terms of radians/sample. We want to convert this to Hz (which is cycle/second). To do so, multiply by Fs/(2*pi).


Time duration of a signal to determine the peak frequency

What window function are you using when you apply your FFT?

If the signal being analyzed is simply truncated on both ends, the window is rectangular.

Lets call the length of time that you generate your signal the analysis window and call the signal you create the "windowed signal".  In this case I'm assuming the window to be rectangular.

One interpretation of the math behind the FFT is that the signal being analyzed is periodic and the complete version of the signal is formed by taking the unaltered input from the analysis window and placing it end to end so it repeats indefinitely as illustrated below

....| Windowed Signal (-1) |  Windowed Signal (0) | Windowed Signal (+1) |...

Depending on where you stop your signal generation, you may introduce discontinuities in the periodic signal that the FFT assumes when you do your analysis.  For example if your signal generation stops at the peak of a sine wave but starts at a 0 crossing, the signal at the point where one windowed signal segment meets the nest will abruptly go from Vpeak to 0.  This discontinuity introduces spurious frequency content in the FFT.  Also if you don't make your analysis window an integer multiple of the period of signal being analyzed, the FFT will assume the signal is periodic in a different period than the signal being analyzed (actually the period assumed by the FFT is always the window length).  this will create frequency lobes that are not harmonically related to the signal you are testing.

As you increase the length of your window, the contribution of the discontinuity and window period diminish, so the longer duration signal will give more accurate results.

The way to improve accuracy of the FFT frequency and/or amplitude resolution is with the use of non-rectangular window functions.  There are a bunch of window functions that are designed for this purpose.  The window function is multiplied by your time domain data (signal).  The effect of the window function is to taper off the amplitude of at the ends of your analysis window (the start and end of the signal you generate) so the edge effects are reduced in the FFT computations.  Different windows have different advantages in frequency or amplitude resolution, so there are trade offs in choosing a window function.

In general, if you want decent resolution in your FFT, you should provide a data length that is least 2x the period of the lowest frequency to be represented and apply a window function to soften the effects of signal truncation.  Note that the longer you make your time window, the greater the frequency resolution of your FFT.  The bins are spaced at 1/T and its multiples where T is the duration of your signal.

There is plenty of information on the internet that discusses the different window functions  used in FFT analysis.  Just Google "FFT window functions".


This problem is related to deconvolution and equalization. You are basically undoing the effect of a filter by another filter, such that the total system has a flat response, i.e., has a unit impulse as its impulse response.

From $$(h\star g)[n]=\delta[n]\tag{1}$$ it follows that $$H(z)G(z)=1\tag{2}$$ must be satisfied. So the solution to the problem is
$$G(z)=\frac{1}{H(z)}\tag{3}$$

as you've already found out.

However, the problem with $(3)$ is that the corresponding filter might not be causal and stable. Note that the poles of $G(z)$ are the zeros of $H(z)$, and the pole locations determine whether $G(z)$ is causal and stable. The inverse filter $G(z)$ is stable if all its poles are inside the unit circle. I.e., all zeros of $H(z)$ must be inside the unit circle for the inverse filter $G(z)$ to be causal and stable. Such a filter $H(z)$ with all its zeros inside the unit circle is called minimum-phase.

As a simple example take $h[n]=\delta[n]+\frac12\delta[n-1]$. We have $H(z)=1+\frac12z^{-1}$, and the zero of $H(z)$ is $z_0=-\frac12$. We can conclude that the inverse filter $(3)$ is causal and stable, because its pole is inside the unit circle:

$$G(z)=\frac{1}{1+\frac12z^{-1}}\tag{4}$$

The inverse $\mathcal{Z}$-transform of $G(z)$ is

$$g[n]=\left(-\frac12\right)^nu[n]\tag{5}$$

It's straightforward to show that the sequences $h[n]$ and $g[n]$ indeed satisfy $(1)$.

In practice, we usually assume that an exact inversion $(3)$ is not possible, so we try to approximately invert a given system. And we usually allow a certain delay, so $(1)$ is changed to $$(h\star g)[n]\approx\delta[n-N]\tag{6}$$ with some appropriately chosen $N>0$.

As a final remark, you have to work on your $\mathcal{Z}$-transform skills. Your attempt at inverting $G(z)$ was not successful.


How to extract features of individual spikes in non-negative signal?

Continuous wavelet transforms can provide information and parameters for sparse, piece-wise regular signals. One example related to yours is present in The Continuous Wavelet Transform in MRS, A. Suvichakorn, C. Lemke, A.  Schuck, J.-P. Antoine, as exemplified in the following picture.


If you were attempting to modify the contents of a DFT before the inverse DFT, say for filtering, then you need to do overlap-add or overlap-save fast convolution instead of just an Inverse DFT.  You will also need to use a longer DFT and inverse DFT.  If you can determine the length of the impulse response of your filter, then the length of the new DFT/IDFT needs to be extended and zero-padded by the length of that impulse response, so that the filter convolution is equivalent to linear, not circular.

If filtering, the noise you hear is likely the result of doing a circular convolution, instead of zero-padding and overlapping.

If, instead you are synthesizing audio, look up phase vocoder methods.


I'm using a cosine-roll-off filter with a roll-off-factor of between 0.2 to 0.5.

Shamelessly copying from wikipedia's Raised-Cosine page, here's a plot of frequency responses of Raised-Cosines filters with different roll-offs:

As you can see, the cutoff frequency for all of these is always half the symbol rate, $\frac12 f_S = \frac12 \frac1{T}$, but the amount of energy outside of that depends on the roll-off. The closer you get to the perfect sinc shape in time ($\beta\rightarrow0$), the sharper that gets.

Notice that this leads to a (double-sided) bandwidth between the cutoff frequencies of $B=2\frac12f_S=f_S$!

However, even for the extreme case of a roll-off factor (which, by the way, is exactly defined to be the difference in occupied bandwidth compared to the minimal bandwidth given by above cutoff frequency) of $\beta=1$, you'd only get $B=2f_S$, so you actually consider the worst case for raised cosine filters – which is a justifiable simplification.

Notice that I'm presuming the book you're citing has a bandwidth definition that covers only a single side!


Calculation of line of sight system gain

You need to multiply the antenna gains, not add them. Specifically, if the free-space loss (attenuation) is $L_{FS}$, the transmitter antenna has gain $G_T$, and the receiver antenna has gain $G_R$, then the total system loss $L$ is $$ L = \frac{L_{FS}}{G_T G_R}. $$ The system gain $G$ is $$G = \frac{1}{L} = G_{FS}G _T G_R, $$ where $G_{FS}$ is the free-space gain.

Of course, if you're doing the calculation in decibels, then the antenna gains are added: $$ G_{dB} = G_{FS,dB} + G_{T,dB} + G_{R,dB}. $$


Least Squares Re Synthesis of Short Time Fourier Transform (STFT)

The least-squares re-synthesis procedure is very similar to the overlap-add (OLA) procedure.

Let w_n be the discrete window vector and y_n(:,k) be the time domain vector computed using the IFFT on a column k of the 2D matrix. Both w_n and y_n(:,k) are the same length.

Then, using Matlab syntax, we compute the point-by-point multiplication with the window:

w_n .* y_n(:,k)

As mentioned in the comments above, the y_n is trimmed to be the same length as the window w_t.  The w_n .* y_n(:,k) is then overlap-added in the same fashion as shown in the code associated with my previous post on the inverse STFT here.  The overlap-added sequence is the numerator of the expression shown in the original question above.

The same overlap-add operation is applied to the squared window w_t.^2.  This is the denominator of the expression shown in the original question above.

Then, the final re-synthesized output is simply the point-by-point division of the numerator by the denominator.  In Matlab syntax,

numerator ./ denominator


How do I construct input to neural network from audio signals?

So the problem is that you want to classify some audio samples but each sample has different sizes. AFAIK there is no special way for doing this with (simple) neural networks.

For general classifiers I am sure there are lots of ways but these are some very simple ones so you can begin playing with the problem:

Divide your signal by frames/segments of equal size and use each frame as if it were a training example. At the end you will use a voting strategy. Since you are using a neural network, you can use the probabilistic outputs of the last layer instead of the hard classes to weight this voting.

Using the average across all frames/segments in your signal. This is a very light-weight approach but can have good results sometimes.

Either way in the "feature extraction" part you will want to use some spectral-based transformation. I recommend either the logarithm of the spectrum or computing MFCC vectors. You can also train your network from raw audio but this is a bit unusual.

You should note that these two methods completely spoil the dynamics of your system because of the granularity. Try playing with the frame sizes or aggregating neighbor frames to increase or decrease the granularity.


The filter doesn't constrain the first sample to be zero. The first sample just has a relatively small value, because it takes the filter some time to react to the input signal. This is because initially the delay elements of the filter contain zeros (i.e. the initial state is zero). This means that for the filter it looks like there has been a long input signal consisting only of zeros, and suddenly the input jumps to a value of around $-20$. This is also why you see the ripples in the output signal. What you see is basically the (negative) step response of the filter. And any frequency selective filter with a decent (frequency domain) performance will have some overshoot in its step response.

Note that the function filter() accepts an initial state as an input argument. So you could compute the filter state for a constant input signal (say about 100 samples or more) equal to the first value of the actual input signal. The filter function will give you the final state after having processed that constant input signal, which you can use as initial state for filtering the actual data. This will greatly reduce any transients in the filter output:

[dum,s] = filter(b,a,x(1)*ones(100,1));    % get initial state
y = filter(b,a,x,s);                       % filter with initial state s

The group delay of the filter is frequency dependent, so it's hard to make general statements about the delay of a specific input signal. What is true is that the average group delay increases with increasing filter order.


What is the effect of reducing NFFT less than the signal length

Reducing the NFFT samples in your DFT will result is less resolution in the spectrum you'll be observing. The frequency resolution is defined as
$$
f_k = \frac{F_s}{N}k\quad \text{with}\quad k=0, 1, \ldots, N-1\tag{1}
$$
In your $N$-point Discrete Fourier Transform (DFT), decreasing $N$ increases the spacing between the $f_k$. Or put in another way, for a fixed sampling frequency $F_s$,

\begin{align}
N_2 &< N_1 \implies \Delta f_{k_2} > \Delta f_{k_1}\\
\text{Where}\quad f_{k_1} &=\frac{F_s}{N_1}k_1\quad\text{and}\quad f_{k_2} =\frac{F_s}{N_2}k_2\\
\implies \Delta f_{k_1} &=\frac{F_s}{N_1}\quad\text{and}\quad \Delta f_{k_2} =\frac{F_s}{N_2}
\end{align}

So say, you have some frequency components (i.e. peaks) you're observing in your spectrum with a high $N$, gradually reducing this value, you'll see that the peaks are getting less and less localized visually. Note that it is only visually, your DFT got your spectrum, you play with $N$ (i.e. NFFT) to visualise more out of your spectrum. Together with $F_s$, it determines the frequency resolution as defined in $(1)$.

EDIT:

To compute the energy, you sum the squared amplitudes spectrum values. Look for instance at the MATLAB code below for  the sine wave
$$
s(t) = 3\sin(2\pi 150t)
$$

Fs = 500;            % Sampling frequency
T = 1/Fs;             % Sampling period
L = 2000;             % Length of signal
t = (0:L-1)*T;        % Time vector
f = Fs*(0:(L/2))/L;

S = 3*sin(2*pi*150*t);

NFFT = 2^nextpow2(L);
NFFT1 = L;
NFFT2 = NFFT/2;
NFFT3 = NFFT/4;

Y = fft(S, NFFT1)/NFFT1;
Y2 = fft(S, NFFT2)/NFFT2;
Y3 = fft(S, NFFT3)/NFFT3;

EY = sum(Y.*conj(Y));
EY2 = sum(Y2.*conj(Y2));
EY3 = sum(Y3.*conj(Y3));

Your FFT have to be normalized and the energy computed as sum of the squared amplitude spectrum. For the sample code above you get:

EY = 4.5
EY2 = 4.49644474953643
EY3 =4.49832143842713

Which is from the amplitude of the sine $3$, giving ($3^2/2$). Also showing the effect of taking different $N$ values for the FFT on resulting energy.


Calculate signal to noise ratio from the Hillbert envelope of a stacked signal

Based on this section of the paper:

Here is my attempt at doing a calculation like that.

The plot above shows three plots:

The average of the dailies (in green) and all the daily data (in black).
The function $\sigma$ calculated from the dailies.
The SNR calculated from $\sigma$ and the dailies.

Scilab code for what I did is below.

// 25624
Tdata = 33;
stack_average = wfir("hp",Tdata,[.3 0],"hm",[0 0]);

Ndays = 30;

ccs = zeros(Ndays,Tdata);

for k=1:Ndays,
    ccs(k,:) =  stack_average + 0.05*rand(1,Tdata,'norm').*window('hm', Tdata);
end;

clf
subplot(311)
plot(ccs','k')
plot(stack_average,'g')
xtitle('$\mbox{All } cc \mbox{ in black and } <cc> \mbox{ in green}$')

sigma_N_t = sqrt((mean(ccs.^2,'r') - mean(ccs,'r').^2)/(Ndays-1));
s_N_t = abs(hilbert(mean(ccs,'r')));
SNR_N_t = s_N_t ./ sigma_N_t;

subplot(312)
plot(sigma_N_t)
xtitle('$\mbox{\rm Level of noise,}\ \sigma$')

subplot(313)
plot(s_N_t, 'r')
xtitle('$\mbox{\rm Signal to noise ratio,}\ SNR$')


Parseval's theorem says that the the following relationship holds
$$
\sum_{n=1}^{N} a[n]\,a^*[n] = \frac{1}{N}\sum_{k=1}^{N} A[k]\,A^*[k]
$$
where $A[k]$ is the discrete Fourier transform of $a[n]$, both assumed to be of length $N$ (no padding). This arises from the fact that the signal energy calculated from the time domain and frequency domain must be equal. See this answer for a little more detail about how this formula comes about.
If you want to calculate a one-sided spectrum, you need to throw away the negative frequencies (which are redundant when the signal is real). However, because half of the energy was contained in these entries, you need to multiply the energy in the remaining bins by 2 - except for the DC and Nyquist. I see that you have recognised correctly that the Nyquist is only present when $N$ is even, and that the first entry is the DC component regardless of $N$ - good. However, you have multipied the amplitudes by 2 (straight after the FFT), not the energy - the negative frequency bins contain half of the energy.
Instead, once you have calculated the discrete Fourier transform using the FFT algorithm (your variables $A$ and $B$), I would first then obtain the energy spectral density of the signal $a[n]$ using
$$
\textrm{ESD}_a[k] = |A[k]|^2 = A[k]\,A^*[k]
$$
and then throw away the negative frequencies from this vector instead. This way you throw away half of the energy correctly, and you can confidently multiply the energy in the remaining bins by $2$.
The following code illustrates this:
%% ODD
a = [1 2 3 4 5];
A = fft(a);
ESD_a = A.*conj(A);
ESD_a_onesided = [ESD_a(1) 2*ESD_a(2:3)];

E_a_timedomain = sum(a.^2)
E_a_twosided = sum(ESD_a)/5
E_a_onesided = sum(ESD_a_onesided)/5

%% EVEN
b = [1 2 3 4 5 6];
B = fft(b);
ESD_b = B.*conj(B);
ESD_b_onesided = [ESD_b(1) 2*ESD_b(2:3) ESD_b(4)];

E_b_timedomain = sum(b.^2)
E_b_twosided = sum(ESD_b)/6
E_b_onesided = sum(ESD_b_onesided)/6

The result is then correctly
E_a_timedomain = 55
E_a_twosided = 55
E_a_onesided = 55

E_b_timedomain = 91
E_b_twosided = 91
E_b_onesided = 91

EDIT----------------------------
Actually, the energy values of $55$ and $91$ obtained are only correct if we are assuming the sampling period of the signal acquisition was $T_s=1$.
The signal energy of a continuous signal $a(t)$ is defined as
$$
E_s = \int_{-\infty}^{+\infty}|a(t)|^2\;dt
$$
and the energy of a sampled version of it is then
$$
E_s = \sum_{n=1}^{N}|a[n]|^2\cdot T_s
$$
and you can see that we need to account for the signal sampling period to get the right energy.
The following verifies Parseval's theorem if we for example had collected the signal at some other sampling period $T_s\neq 1$
Ts = 0.05; % the sampling period of acquisition
Fs = 1/Ts; % the sampling frequency of the acquisition

%% ODD
a = [1 2 3 4 5];
N = 5

A = fft(a)*Ts;
ESD_a = A.*conj(A);
ESD_a_onesided = [ESD_a(1) 2*ESD_a(2:3)];

E_a_timedomain = sum(a.*conj(a))*Ts
E_a_twosided = sum(ESD_a)*Fs/N
E_a_onesided = sum(ESD_a_onesided)*Fs/N

%% EVEN
b = [1 2 3 4 5 6];
N = 6

B = fft(b)*Ts;
ESD_b = B.*conj(B);
ESD_b_onesided = [ESD_b(1) 2*ESD_b(2:3) ESD_b(4)];

E_b_timedomain = sum(b.^2)*Ts
E_b_twosided = sum(ESD_b)*Fs/N
E_b_onesided = sum(ESD_b_onesided)*Fs/N

with the output
E_a_timedomain = 2.75 [signal^2 sec]
E_a_twosided   = 2.75 [signal^2 sec]
E_a_onesided   = 2.75 [signal^2 sec]

E_b_timedomain = 4.55 [signal^2 sec]
E_b_twosided   = 4.55 [signal^2 sec]
E_b_onesided   = 4.55 [signal^2 sec]


The main trick is knowing how typical lumped circuit element impedances are represented in the $s$-domain. Recall the current-voltage relationship for a capacitor:

$$
i(t) = C \frac{dv(t)}{dt}
$$

Transforming this to the Laplace domain yields:

$$
I(s) = CsV(s)
$$

Note that voltage and current still have a proportional relationship in this domain, just like the time domain. This allows us to write a (frequency-dependent) $s$-domain impedance for the capacitor. This impedance can be used with Ohm's law and techniques like nodal analysis to calculate what a circuit's Laplace-domain output would be for a given Laplace-domain input signal (e.g. an impulse, whose Laplace transform is unity).

For a capacitor, this impedance is:

$$
Z(s) = \frac{V(s)}{I(s)} = \frac{1}{Cs}
$$

This gives capacitors their known highpass characteristic; for $s=j\omega$, as $\omega$ gets larger, the impedance gets smaller. For the limiting case of $\omega=0$, $Z(0) = \infty$; a capacitor passes no DC current.

We can follow the same procedure for an inductor, starting with its current-voltage relationship:

$$
v(t) = L\frac{di(t)}{dt}
$$

Apply the Laplace transform:

$$
V(s) = LsI(s)
$$

and again, we can calculate the inductor's $s$-domain impedance:

$$
Z(s) = \frac{V(s)}{I(s)} = Ls
$$

As the dual to capacitors, inductors have a lowpass characteristic. For small frequencies ($s \to 0$), $Z(s) \to 0$ as well.

Given these impedances, it becomes very straightforward to apply basic circuit analysis techniques to determine the $s$-domain transfer function of a network containing resistors, capacitors, and inductors.


why LTE doesn't use DC subcarrier in the downlink?

The DC subcarrier is normally omitted because in a zero-IF receiver it overlaps with DC offsets in the hardware which are very difficult to remove.


Measuring the average frequency of signals

The weighted average "(sum over (fourier coefficient index * amplitude of that coefficient)) / number of fourier coefficients" does not seem very good as it will always give more weightage to the higher frequencies( because fourier coefficient index for higher frequencies will be higher) and the frequencies which you are looking for might be having low index.

Better method could be to see which fourier coefficients have maximum values for each bang of hammer and corresponding to that you can find the major frequencies(band of frequency) for each bang.

hope this helps.


Here is one way to think about this:

Start with the rectangular pulse of 10ms length. DFT is the sinc() function with a "period" of 100Hz
Now repeat periodically with 100ms intervals. Mathematically that's a convolution with a 100ms pulse train.
The DFT of the pulse train is also a pulse train in the frequency domain with 10 Hz spacing
Convolution with in the time domain is multiplication in the frequency domain. Hence we need to multiply the 100 Hz sinc() with a 10 Hz pulse train. That's the basic shape you have: a line spectrum with 10 Hz line spacing with the overall contour of a 100Hz sinc()
Sample at 1kHz. The sinc() function decays fairly slowly with frequency: around the Nyquist frequency you are only at -25 dB or so. Hence you will get non-trivial aliasing. However, since your sample rate is an integer multiple of the line spacing, the aliasing products will fall on top of the existing lines and will be hard to see. It'll just be a small deviation from the sinc() contour.

The "fuzz" between the lines appears to be noise. Your time domain wave form looks noisy and not particularly clean. It could also be caused by the way you calculate the spectrum (windowing, periodicity misalignment, FFT size, etc.)


Multiple choice on sampling and aliasing

Question 1
Perhaps the easiest way to understand this is to draw $X(j\omega)$:

With the corresponding aliases after sampling in blue:

So long as the blue rectangle do no overlap you would not have aliasing. At this point it should be more clear that you need
$$
\omega_s \gt \omega_0 \\
T_s \lt \frac{2\pi}{\omega_0} \\
$$
Question 2
Again, plotting $G(j\omega)$ (red), together with $75X(j\omega)$ and an alias image after sampling by $\omega_s$ (dashed blue) can help visualize the situation:

Indeed, because $\omega_{\max} > \frac{1}{2}\omega_s$ you have some aliasing. However the question is not whether you can avoid aliasing, but rather to find the value $\omega_0$ below which the aliasing does not affect the spectrum.
In the above graph you may notice that the red and dashed blue curves overlap up to $50\pi$. So $\omega_0 = 50\pi$.


What you actually want to do is the following: you need a zero at $\omega=\pi/4$, but you also need to place a zero at $\omega=-\pi/4$ in order to get real-valued filter coefficients. So your transfer function becomes

$$\begin{align}H(z)&=\frac{(z-e^{j\pi/4})(z-e^{-j\pi /4})}{z^2}\\&=\frac{z^2-2z\cos(\pi/4)+1}{z^2}\\&=1-\sqrt{2}z^{-1}+z^{-2}\tag{1}\end{align}$$

The corresponding difference equation is

$$y[n]=x[n]-\sqrt{2}x[n-1]+x[n-2]\tag{2}$$


Laplace Transform of $f(t+a), a>0$ where $f(t)$ is not periodic

Let $s = \sigma + j\omega$, the inverse Laplace transform of $f(t+a)$ is given by
$$f(t+a) = \frac{1}{2\pi j} \int_{\sigma-j\infty}^{\sigma+j\infty} F(s)e^{s(t+a)} \mathrm{d}s = \frac{1}{2\pi j} \int_{\sigma-j\infty}^{\sigma+j\infty} F(s)e^{sa}e^{st} \mathrm{d}s.$$

Hence the bilateral Laplace transform of $f(t+a)$ is $F(s)e^{sa}$ where $F(s)$ is the Laplace transform of $f(t)$. For the unilateral case, see Matt L. answer.

This is sometimes called the shifting theorem. See Theorem 12.16 here.


Calculation of inverse impulse response in the frequency domain

Your result is correct, even if it seems a bit counterintuitive. Note that you don't compute the Fourier transform but the discrete Fourier transform (DFT). The system you want to invert has the impulse response
$$g[n]=\frac12\delta[n-19]\tag{1}$$
Obviously, the inverse system must advance the incoming signal by $19$ samples:
$$g_{inv}[n]=2\delta[n+19]\tag{2}$$
Clearly, that system is non-causal. Taking into account the implicit periodicity of the DFT (and IDFT) with length $N=100$, the signal $g_{inv}[n]$ cannot be distinguished from a signal $g_{k,inv}[n]=g_{inv}[n+kN]$, $k\in\mathbb{Z}$. Your result is just $g_{-1,inv}[n]$:
$$g_{-1,inv}=2\delta[n+19-N]=2\delta[n-81]\tag{3}$$


Why then is the pulse shaping filter not used to filter the sent signal in order to draw the BER curve?

When drawing a BER curve, we're considering the bit or symbol energy.
That energy is the convolution of transmit pulse shape with receive filter evaluated at the symbol instance.
So, that filter is taken into account. We just usually define it such that the pulse energy is 1, and then you don't see it in your equations.


Transfer function of a PLL Loop Filter that can support a linearly increasing (chirping) frequency

To track a frequency ramp with a Phase lock loop, with zero steady state error requires a type 3 PLL Loop; which means three integrations (DC Poles) in the open loop gain (your NCO would be one of the integrators and your loop filter needs to provide the other two). Stabilizing such a system becomes more challenging but here is one reference paper detailing an approach to the design:

http://vbn.aau.dk/files/72548181/Final_PLL.pdf

This picture helps clarify tracking and loop order; given a phase lock loop, phase is the parameter being tracked. A frequency ramp as you are dealing with is an accelerating phase, hence the requirement for a type 3 loop.  The detailed loop filter design would depend on your overall system and related gain coefficients for your specific design, but I hope the above reference offers some guidance to help you.

Also note that many people often confuse "Type" with "Order", so to be very clear a Type Three Loop is not the same as a Third Order Loop. Specifically a Type Three Loop must be Third Order (or higher), but a Third Order Loop need not be Type 3. Type specifically is number of integrators in the loop, which are poles at z=1 in the open loop gain, while order is the total number of poles, regardless of location (specifically the order of the denominator polynomial of the open loop gain, which ends up being the order of the closed loop gain given proper transfer functions).


Output of an LTI system given its transfer function and input

You cannot solve this problem using the Laplace transform. The reason is that the Laplace transform of the input signal doesn't exist. You could use the Fourier transform, but in this case there's an even simpler way to determine the output signal. You need to know one important property of linear time-invariant (LTI) systems: their response to a sinusoidal input is a sinusoidal signal with the same frequency, but with its amplitude and phase altered according to the system's frequency response evaluated at the input frequency.

So for an input signal

$$x(t)=A\sin(\omega_0t+\phi)\tag{1}$$

the output is given by

$$y(t)=A\big|H(j\omega_0)\big|\sin\left(\omega_0t+\phi+\arg\big\{H(j\omega_0)\big\}\right)\tag{2}$$

where

$$H(j\omega)=\big|H(j\omega)\big|e^{j\arg\{H(j\omega_0)\}}\tag{3}$$

is the system's frequency response.


Overlap/add time-domain audio frames: How does normalization/scaling work with overlap greater than 50%?

Let's assume continuous time (rather than discrete time).

If you do not process the windowed data at all, you would like the output (the sum of the windowed frames) to be equal to the original signal. Allowing scaling of the output by a constant scaling term, this is only possible if the sum of all of the time-shifted window functions is constant over time.

A typical window function (including Hann) is a sum of harmonic cosine terms. Let's assume you only use such windows. At the end of the window, the cosine terms have derivative discontinuities that can only be "cured" by stitching to the end of the window the start of another window which has identical values of all derivatives. This, and that you probably want the windows to be distributed evenly, requires that the hopsize is $\frac{1}{N}$ window lengths, with $N$ integer. The truncated cosine terms in each of the first $N$ windows are extended to complete cosines by concatenation of every $N\text{th}$ window.

The sum $f(T)$ of the windows distributed this way is:

$$f(T) = \sum_{n = 0}^{N - 1}\sum_{h = 0}^{H - 1} a_h\cos\left(2\pi h\left(T-\frac{n}{N}\right)\right),$$

where $h$ is the harmonic number, $H$ is the number of harmonics in the window function, $a_h$ are the window function specific coefficients, $T$ is time expressed so that a time difference of $1$ means one window length. $h$ and $n$ are the running harmonic and hop numbers.

Let's try window densities $N = 1\dots3$ with a window that has two harmonics ($H=2$). First $N=1:$

$$f(T) = a_0 + a_1\cos(2\pi T)$$

then $N=2:$

$$f(T) = a_0 + a_1\cos(2\pi T) + a_0 + a_1\cos\left(2\pi \left(T - \frac{1}{2}\right)\right) = 2a_0$$

and $N=3:$

$$f(T) = a_0 + a_1\cos(2\pi T) + a_0 + a_1\cos\left(2\pi \left(T - \frac{1}{3}\right)\right) + a_0 + a_1\cos\left(2\pi \left(T - \frac{2}{3}\right)\right) = 3a_0$$

This could be for example the Hann window with $a_0 = \frac{1}{2}$ and $a_1 = \frac{1}{2}:$

Figure 1. Hann windows sum to constant with a hopsize of 1/2 or 1/3 window lengths, but not with a hopsize of 1 window length.

Although it is not as evident as with the 1/2 hopsize, the $1\text{st}$ harmonic cosines cancel each other completely with the 1/3 hopsize. The cancellation can be illustrated by representing the phases and amplitudes of the harmonic terms as complex numbers:

Figure 2. Complex number representation of the amplitudes and phases of the 1st harmonic terms in shifted concatenated Hann windows, with a hopsize of $1/3$ window lengths. The rotational symmetry of the pattern means that the first harmonics cancel, because the sum of the three complex numbers must be zero.

For window functions with $H$ harmonic cosine terms, if we tabulate for different window densities $N$ the factors appearing in front of each $h\text{th}$ harmonic cosine term in $f(T),$ a pattern emerges:

$$\begin{array}{r|lllllllllll}&\rlap{h}\\
N&0&1&2&3&4&5&6&7&8&9\\
\hline
1&1&1&1&1&1&1&1&1&1&1\\
2&2&0&2&0&2&0&2&0&2&0\\
3&3&0&0&3&0&0&3&0&0&3\\
4&4&0&0&0&4&0&0&0&4&0\\
5&5&0&0&0&0&5&0&0&0&0\\
6&6&0&0&0&0&0&6&0&0&0\\
7&7&0&0&0&0&0&0&7&0&0\\
8&8&0&0&0&0&0&0&0&8&0\\
9&9&0&0&0&0&0&0&0&0&9\end{array}$$

The factors in the 2-term window example appear in the two first columns and three first rows.

For $H \le N,$ all harmonics except for the $0\text{th}$ harmonic cancel completely. This means for example that the 4-cosine-term ($H = 4$) Blackman–Harris window could be used with a 1/4 hopsize or greater ($N \ge 4$), requiring only an output scaling by $\frac{1}{4a_0} = \frac{1}{4\times0.35876}.$


noise cancellation of non-repetitive noise

Since I cannot comment, here comes a short answer.

As far as I know, there are microphones outside the earphone. The microphones record the outside noise and generate sounds with 180 phase shift to cancel the noises. The noises are canceled actively.


Why eigen values and poles of a system are equivalent?

Let's consider a discrete-time state space model (the derivation for a coninuous-time system is completely analogous):

$$\begin{align}\mathbf{q}[n+1]&=\mathbf{Aq}[n]+\mathbf{b}x[n]\\
y[n]&=\mathbf{c}^T\mathbf{q}[n]+dx[n]\tag{1}
\end{align}$$

where $x[n]$ is the input, $y[n]$ is the output, and $\mathbf{q}[n]$ is the state vector. Taking the $\mathcal{Z}$-transform of (1) gives

$$\begin{align}z\mathbf{Q}(z)&=\mathbf{AQ}(z)+\mathbf{b}X(z)\\
Y(z)&=\mathbf{c}^T\mathbf{Q}(z)+dX(z)\tag{2}
\end{align}$$

From the top equation of (2) we can express $\mathbf{Q}(z)$ in terms of $\mathbf{A}$, $\mathbf{b}$, and $X(z)$:

$$\mathbf{Q}(z)=(z\mathbf{I}-\mathbf{A})^{-1}\mathbf{b}X(z)\tag{3}$$

where $\mathbf{I}$ is the identity matrix of appropriate dimension. Plugging (3) into the second equation in (2) gives

$$Y(z)=\mathbf{c}^T(z\mathbf{I}-\mathbf{A})^{-1}\mathbf{b}X(z)+dX(z)\tag{4}$$

from which follows the expression for the transfer function

$$H(z)=\frac{Y(z)}{X(z)}=\mathbf{c}^T(z\mathbf{I}-\mathbf{A})^{-1}\mathbf{b}+d\tag{5}$$

The matrix $(z\mathbf{I}-\mathbf{A})^{-1}$ can be written as

$$(z\mathbf{I}-\mathbf{A})^{-1}=\frac{1}{\text{det}(z\mathbf{I}-\mathbf{A})}\mathbf{B}\tag{6}$$

where $\mathbf{B}$ is the adjugate matrix of $(z\mathbf{I}-\mathbf{A})$ whose elements are polynomials in $z$.

Plugging (6) into the expression for the transfer function (5) shows that the denominator of $H(z)$ is given by $\text{det}(z\mathbf{I}-\mathbf{A})$, which is the characteristic polynomial of $\mathbf{A}$ the roots of which are the eigenvalues of $\mathbf{A}$. Since the roots of the denominator of $H(z)$ are the system's poles, the eigenvalues of $\mathbf{A}$ and the poles of $H(z)$ are equivalent.


It seems like you don't want to make any prior assumptions about the exact shape of the signal, other than assuming that it is ``smooth''. This is best attacked using a non-parametric function fitting/regression approach. Use either
a) spline fitting or b) kernel smoothing.

Two reasons why this is a good idea:

These methods are optimal in the sense that no other non-parametric filtering technique (moving average, etc.) can give you a better worst case error decay rate. This answers your main question: there is no better method if you care about minimax error.
Spline fitting and kernel regression codes are very easily available.

A quick Google search gives me the following:
In Octave interp1 has a 'cubic' option:
http://www.gnu.org/software/octave/doc/interpreter/One_002ddimensional-Interpolation.html

Open source code for kernel smoothing:
http://code.google.com/p/pmtk3/source/browse/trunk/matlabTools/stats/kernelSmoothingRegression.m?r=2645&spec=svn2776


Extracting original signal from overlapping area

I'd calculate the horizontal point spread function for each color channel separately and use the spectral sensitivity ratios of the channels for frequency calibration; the resulting sensitivity-weighted spectra should obey the same ratios.

Deconvolution can be implemented as frequency domain division of each diffracted row by the corresponding direct transmission row. The division results can be averaged over all rows, excluding from statistics any noisy bins that came from division by a value too close to zero. Hopefully the rows differ enough to cover all frequency bins. Combining the color channels at the end further reduces noise.

For this to work properly, the image should be corrected for any distortion and rotated so that the point spread function due to diffraction is fully horizonal. Chromatic aberrations will be hard to fix.

A completely different approach would be kind of a matching pursuit, which should work quite fast thanks to the spiky spectra of fluorescent phosphors. There you'd try to pattern-recognize what position in the diffracted image is most similar to an intensity-scaled copy of the direct transmission image (using a color appropriate to that position, according to a calibration). Then you'd subtract the match from that position in the diffracted image and simultaneously accumulate an estimate of the spectrum. You'd repeat the process until you'd get good enough an estimate, or would no longer be able to find a good match. It should be possible to do image geometry correction in the pattern recognition part of this approach, which would also be able to avoid problems due to chromatic aberration.


What does it mean to plot a PAM signal over a number of symbols?

$10$ symbols means that you choose $N=10$, i.e., you plot realizations of

$$u_c(t)=\sum_{n=0}^{9}b_c[n]p(t-n)$$

and

$$u_s(t)=\sum_{n=0}^{9}b_s[n]p(t-n)$$

in the range $t\in[0,10]$.


Comparing spectrogram similarity over time

For spectrogram similarity, I suggest using a metric from image processing known as SSIM - Structural Similarity Index Metric. It is a quantity between $[0,1]$ which denotes the similarity between two images.
The paper is called Image Quality Assessment: From Error Visibility to Structural Similarity. The MATLAB function ssimval() will give you the metric.


How do I eliminate line noise from USB microphone

You already know the noise frequencies, if you know the phases, you can cancel the noise by adding mods with opposite phase, like adaptive headphones do.


How to measure the quality (or difference in quality) of different images?

It all depends on what you mean by "quality".  If it is visual quality as perceived by people, then the best thing to do is to ask multiple people to compare pairs of images to see which ones they like better.

Since you have mentioned a vision application, then it must have some performance metrics for whatever it does, which would be its measure of quality.

If by quality you mean how blurry or noisy the image is, then you might compute the Fourier transform to see whether the image contains high frequency components to see how blurry it is, or compute its entropy to see how noisy it is.


Recovering time domain signal from sum of independent and identically distributed signals

The answer was given in a comment in a cross-post. I wasn't sure if it was a computational or conceptual issue I was having. Turns out it just an issue with the way FFT expects the input to be ordered and how the output is sorted.


Discrete Convolution of two piecewise sequences having this specific form

Given your Eq-3
$$
(x * h)[n]=\sum_{k=\max (a, n-d)}^{\min (b, n-c)} x[k] h[n-k] . \tag{1}
$$
where $x[n]$ and $h[n]$ are as in Eq-1, then you will have
$$
\sum_{k=\max (a, n-d)}^{\min (b, n-c)} \alpha \beta = \alpha \beta (\min (b, n-c)- \max (a, n-d) + 1). \tag{2}
$$
Where you should evaluate the closed form expression Eq.2 for the non-zero range of output given by the range of $n$ :
$a+c \leq n \leq b+d $


GMSK Big error in demodulation MatLAB Simulink

That’s is because of the delay in the GMSK demodulator. MATLab GMSK demodulator uses Virtebi algorithm. You should delay the Rx BER in order to get back the correct data alignment.


Doubt on the relationship between the sensitivity and the complementary sensitivity function

You already have the answer $S+T = 1$, the relation must hold true at all frequencies.

Therefore if your sensitivity function $S$ has a high-pass characterisitic, the complementary sensitivity function $T$ will have low-pass characteristic.

As for bandwidth, your curves make sense. The red sensitivity function has a larger bandwidth than the green sensitivity function. Therefore, the red complementary sensivitity function has a smaller bandwidth than the green complementary sensitivity function.


Actually you have quite done it, but lets indicate.

When the frequencies of both sinusoidal signals $x_1[n]$ and $x_2[n]$ are the same; i.e., $k_1 = k_2 = m$, then the impulses in the corresponding N-point DFT sequences $X_1[k]$ and $X_2[k]$ will occur at the same index $k=m$ and $k=-m$ (or $N-m$) with a weight of $N/2$.

Hence multiplication of $X_1[k]$ with $X_2[k]$  will be $N/2 * N/2 + N/2 * N/2 = N^2/2$, as there are two impulses in the DFT of a cosine signal...

Finally dividing by $N$ yields the result as $N/2$ when the frequencies are the same and $0$ when not.


Impulse Response of a Moving Average System

Taking Equation $2.74$ from Example $3$ and setting $M_1=0$ gives:
$$ h[n] = \frac{1}{M_2+1} \sum_{k=0}^{M_2} \delta[n-k] $$
Now lets take a closer look at this: $h[n]$ is non-zero only for certain values of $n$. Its a good a starting point as any, so lets look at $n=0$. We have:
$$ h[0] = \frac{1}{M_2+1} \bigg( \delta[0] + \delta[-1] + ... + \delta[-M_2] \bigg) = \frac{1}{M_2+1}$$
For $n=1$, we have:
$$ h[1] = \frac{1}{M_2+1} \bigg( \delta[1] + \delta[0] + ... + \delta[1-M_2] \bigg) = \frac{1}{M_2+1}$$
For $n=2$, we have:
$$ h[2] = \frac{1}{M_2+1} \bigg( \delta[-2] + \delta[-1] + \delta[0]+ ... + \delta[2-M_2] \bigg) = \frac{1}{M_2+1}$$
And this pattern continues, the only non-zero value, $\delta[0]$, percolates through until $n=M_2$, then we have:
$$h[M_2] = \frac{1}{M_2+1} \bigg( \delta[-M_2] + \delta[-M_2+1] + ... + \delta[0] \bigg) = \frac{1}{M_2+1}$$
And for $n = M_2 + 1$:
$$ h[M_2+1] = \frac{1}{M_2+1} \bigg( \delta[-M_2-1] + \delta[-M_2] + ... + \delta[-1] \bigg) = 0$$
And you can see, for all $n > M_2$, $h[n]=0$. So, altogether we can write $h[n]$ as equal to $\frac{1}{M_2+1}$ for $n \in [0, M_2]$ and zero otherwise.
All Equation $2.90$ does is treats the summation of impulses (from Equation $2.74$) as unit step functions, a very useful trick! Since $h[n]$ starts at $n=0$ and ends at $n=M_2$ we do: $u[n]-u[n-M_2-1]$ to get a length $M_2$ box car, then scale it by $\frac{1}{M_2+1}$ to get the right amplitude. Hope this helps!


How to determine the phase of a time-sampled periodic signal?

You could try to Fourier transform,

$F(k) = \int_{-\infty}^{\infty} x(t) e^{-2\pi i kt} dt$

which can be expressed in terms of magnitude and phase in the following way

$F(k) = |F(k)|e^{i\Phi(k)} = a(k) + ib(k)$

with magnitude $|F(k)|$ and phase $\Phi(k) = \tan^{-1}(\frac{b}{a})$ for each frequency component $k$. Since it looks like you have a dominant $k$ representing your sine wave you can find the phase of that wave pretty easily. There are probably more obvious ways than this though.


What is the relation between the carrier frequency and the attenuation of electromagnetic waves?

60GHz is associated with absorption by Oxygen.

Rain and Snow reduce range.

Radio waves can penetrate fresh water.  The Great Lakes was a place where submarine crews trained during WW2 and could send and receive underwater. Salt water is another situation entirely. Frequencies need to be very low.  There is often a good natural wave guide just above the ocean surface.

Foliage can attenuate.

The situation in the optical and IR regions are complicated. There is a program called MODTRAN that models absorption in detail.

For 60GHz, attenuation can actually be a virtue, not a problem. You can pack more receivers and transmitters closer together, at a given complexity.


Why do I see a frequency shift of the peak between the S11 and impedance response?

[This is more of an electronics related question and perhaps should be migrated to Electronics.SE?]
This can occur with a series parasitic inductance or capacitance in the resonator; the case with a series inductance is demonstrated below.
The Smith Chart visual provides immediate insight of a shunt L-C resonator (tank circuit) to ground. In the first plot where we have no parasitic inductance in series with an exaggerated parasitic resistance in the resonator (reducing the input impedance at resonance when the reactive inductance cancels the reactive capacitance). At resonance $X_{Lp}$ cancels $X_{Cp}$ and we would be at the large dot shown (with just the shunt resistance due to finite Q remaining, in parallel with the 50 ohm load and thus reducing the input impedance along the resistance axis). As we move off resonance lower in frequency the tank circuit is inductive and ultimately a short to ground at DC, following the shunt inductance trajectory as shown by the path marked $L_p$ indicating a parallel inductor to ground, as we sweep the test frequency toward DC terminating on the left origin in the Smith Chart consistent with a short circuit. As we go higher in frequency the tank circuit will be capacitive and ultimately a short to ground at $f=\infty$, following the shunt capacitor trajectory as shown by the path marked $C_p$ as we sweep the frequency higher.   In this case for all other frequencies other than at resonance the net impedance of the tank would smaller with a maximum tank impedance at resonance, and the resonance peak would match the S11 null.

Observe if the tank circuit had a series inductance effectively prior to the shunt L and C to ground. This inductor in series would follow the trajectory marked by $L_s$ (and would ultimately end at the right origin where $Z=\infty$, terminating at the large dot shown at resonance with $X_{Lp} = -X_{Cp}$. As we move off resonance upward in frequency, the resonator adds a shunt capacitance that in this case brings us closer to the center origin (where Z=50 ohms, perfect match), providing a lower S11 as it maximizes the net impedance of the resonator's connection to ground; essentially the resonator is tuning out the series inductance as it becomes capacitive. Not detailed here is that series inductance would also increase as frequency increases but we can see how if the capacitance relative to the inductance is significantly larger that the trajectory would move down faster than up and to the right resulting in the minimum as the OP is seeing. Further the inductance itself would have a point where it resonates out and becomes capacitive, so we could also easily envision the possible case of a net circular path of the series inductance and shunt capacitance of the resonator passing right through the center of the Smith chart (or closer to it).

Below is a basic simulation of this compliments of https://www.will-kelsey.com/smith_chart/
At resonance the Smith chart would be at DP5 (where the red line intersects the black line to the upper left of center. As we sweep frequency, the trajectory follows the red line where we see a minimum impedance is reached which is not on the resonance of the shunt L and C.


Hardware requirements for ambient active noise cancallation system using adaptive LMS

I would recommend checking out "Design of active noise control system with the TMS320 family" by Sen Kuo [[report link]].

TMS320 chips are specifically designed to perform various ANC tasks. In the report, they even provide code for programming LMS and FxLMS.

In a lot of the literature on real-time ANC, the control system is implemented using these chip sets. So in my opinion, these chips would be a good place to start. Alternatively, I'm sure this system could be implemented with a variety of microcontrollers, but I am not aware of any other hardware specifically designed for ANC.

Side note: I would consider looking into the leaky FxLMS algorithm as well. It's no more complex than regular FxLMS, but provides more stability to the adaptive algorithm (although reduces the steady state performance). The leakage term is typically necessary in real-time practical control systems.


What is the role of Freqency in Digital Images

You can take images as 2D discrete signals. The "time" in 1D signals is actually two spatial dimensions in images (2D signals).

You can measure "frequency" as well - imagine a line of white pixels with regular spacing. The spacing represents period $p$, and frequency is given by $1/p$. Hence the maximum frequency the discrete signal can contain is limited by the sampling rate.

When creating a digital image, sampling can cause new frequencies to appear due to aliasing. See Nyquist frequency.

Fourier and wavelet transforms have discrete variants to work with discrete signals - very simply said: integrals are replaced by sums.

Wavelets are more efficient in image compression, because they can deal with localized redundancies. Simply look of DFT of a grey square on black background and then on DWT of that square.

DWT will exploit transitions of the square shape in high fequency layer, and other layers will be almost clear. The sharp 1px transition would be best described with Haar wavelet without any remainder.

On the other hand, DFT or DCT will result in many cosine functions which summed together form the edge without visible ringing artifacts (Gibbs effect).

Wavelet have its artifacts too, but these are better localized, hence does not corrupt rest of the image (or rest of the macroblock in case of JPEG/DCT compression).


If you can do it that way, it makes it straight forward to implement the Kalman filter. As an example - consider a constant velocity model in Cartesian coordinates, but the only thing you are measuring is the position (not velocity) and you are measuring the position in a Cartesian coordinate system. Then your $\mathbf{H}$ matrix will only pick off the position value as you've described.

Now consider a Sonar or Radar system - which measure range and bearing rather than x, y, z. Now the measurements take place in a Polar/Spherical coordinate system but it is easier to keep the model for your targets as having constant velocities in the cartesian coordinate system. Now your state propagation matrix is the same as in the previous case, but now your $\mathbf{H}$ matrix needs to convert your Cartesion state coordinates to Polar/Speherical coordinates for your measurements. e,g, $R =\sqrt{x^2+y^2+z^2}$. So now your $\mathbf{H}$ matrix is not simply extracting one of the state variables. In fact your $\mathbf{H}$ is not even a constant value - it varies with the values of the state vector. Thus the transformation in the measurement matrix $\mathbf{H}$ is non-linear, so the Kalman filter is not even optimal anymore. In this case you now get into the Extended Kalman Filter to try and handle the non-linearity in the transformation.

This was meant to give you one specific example. Also, note that there are other ways of handling the Cartesian to Spherical conversion issue.


Finding the frequencies given in a wav file

Polyphonic pitch estimation is still a difficult problem.  Search on the MIREX forum for current research papers.  Simplistic spectrum based approaches usually require component pitches with extremely simple timbres (nearly pure unmodulated sine waves, etc.) rarely found in real-life polyphonic music.  More typical polyphonic music can contain a dense set of overlapping overtones plus other transient noises.  More complicated timbres creating polyphonic chords may require a statistical pattern matching based approach (against a database of potential instrument sounds), and/or machine learning algorithms, in addition to any signal processing for spectrum and envelope information.


Which step response matches the system transfer function

The final value of the step response is the DC gain of the closed-loop transfer function, which is generally different from the open-loop DC gain.
Assuming unity gain feedback, the feed-forward transfer function $G(s)$ equals the open-loop transfer function, and the closed-loop transfer function is given by
$$C(s)=\frac{G(s)}{1+G(s)}\tag{1}$$
The final value of the closed-loop step response equals $C(0)$. With the given value of $G(0)$, which I would interpret as $-3$dB, the final value of the closed-loop step response should be around $\sqrt{2}-1\approx 0.4$. The step response in Figure $1$ seems like a good match.


Image Processing: How to detect a quadrialteral signboard in image?

Both your images contain many lines that have nothing to do with the sign you're looking for. And some of those lines are longer/have higher contrast than the lines you actually want, so I think detecting the edge lines (e.g. using a hough transform or by summing up contrasts horizontally/vertically) will not work.

But: The sign you're looking for has other characteristics that should be easier to detect:

There sign background has (almost) constant brightness
It takes up a relatively large area of the image
It's near the center of the image

So you're looking for a large connected area with low contrast. I've hacked a proof-of-concept algorithm in Mathematica. (I'm not an OpenCV expert, but I'll mention the respective OpenCV function when I know them.)

First, I use gaussian derivative filters to detect the gradient magnitude at each pixel. The gaussian derivative filter has a wide aperture (11x11 pixels in this case), so it is very noise-insensitive. I then normalize the gradient image to mean=1, so I can use the same thresholds for both samples.

src = Import["http://www.freeimagehosting.net/uploads/720da20080.jpg"];
pixels = ImageData[ColorConvert[src, "Grayscale"]];
gradient = Sqrt[GaussianFilter[pixels, 5, {1, 0}]^2 + GaussianFilter[pixels, 5, {0, 1}]^2];
gradient = gradient/Mean[Flatten[gradient]];

OpenCV implementation: You can use sepFilter2D for the actual filtering, but apparently, you'll have to calculate the filter kernel values yourself.

The result looks like this:

In this image, the sign background is dark and the sign borders are bright. So I can binarize this image and look for dark connected components.

binaryBorders = Binarize[Image[gradient], 0.2];
sign = DeleteBorderComponents@ColorNegate[binaryBorders];
largestComponent = SortBy[ComponentMeasurements[sign, {"Area", "ConvexVertices"}][[All, 2]], First][[-1, 2]];

OpenCV implementation: Thresholding should be straightforward, but I think OpenCV doesn't contain connected component analysis - you can either use flood fill or cvBlobsLib for that.

Now, just find the largest blob near the center of the image and find the convex hull (I've simply used the largest blob that's not connected to the background, but that might not be enough for every image).

Results:


How to eliminate negative frequencies from IQ signal

I managed to do this now by having a 3 step process after converting the I/Q data to complex numbers. This is far from me understanding what i'm doing, but it's eliminating negative frequencies.

I shift the the frequencies by multiplying each i/q sample with the complex output of an oscillator at frequency PI/2. The oscillator outputs a complex number with Real being sin(pi*0.5*t), and Imaginary being cos() of the same values. Then I multiply that complex number with each sample (t advances +1 with each sample).
Next step is to high-pass filter at a quarter sample rate. So if our sample rate is 192000, my high pass cutoff is at 48000.
Lastly, we shift the whole thing back by doing the same multiplication as in 1), except with PI*1.5.

The result is I/Q data that's got everything in the negative spectrum filtered out.

I read about Hilbert transforms etc etc, but got lost halfway through that. Maybe someone has an explanation of what i did or why this works? :)


How to use values from the output of a block in GNURadio

I want to use the output values of the top_block

A top_block has no output; a signal processing block in a flowgraph has.

So I want the flowdiagram to take one path if the output of my processing block is 0 and take a different path if it is non-zero.

Ah, you kind of want to "demux" something. So write a simple demux block that has two inputs and multiple outputs; copy from the first input to any of the outputs, depending on your second input.

I have implemented an if statement within the top block which disconnects and connect parts of the diagram accordingly.

Don't do that. Flow graph reconfiguration is asynchronous to signal processing in GNU Radio, i.e. this won't work, and especially not on a per-sample basis.


Microphone Array null steering unexpected results after IFFT

I don't understand what you did but I think you didn't understand the idea correctly.
Lets say the weighting vector at a certain frequency is given by $\mathbf{w}^H$ where ${}^H$ is the complex conjugate. That is the elements of $\mathbf{w}^H$ are the weights you apply to each sensor without null steering. And lets say that the array response to the interference is $\mathbf{v}_I$ (which is basically the conjugate of the array weights if you were to point the beamformer to the interference). The simple null steering solves the following optimisation problem:
$$   \underset{\mathbf{w}_{d}}{min}\:\left\Vert \mathbf{w}_{d}-\mathbf{w}\right\Vert ^{2}\\
  subject\; to\\
  \qquad\mathbf{w}_{d}^{H}\mathbf{v}_{I}=0
$$

Using Lagrange multipliers it can be shown that
$$ \mathbf{w}_{d}^{H}=\mathbf{w}^{H}\cdot\left(\mathbf{I}_{N}-\frac{1}{N}\mathbf{v}_{I}\mathbf{v}_{I}^{H}\right)$$
where $\mathbf{I}_{N}$ is the identity matrix and $N$ is the number of sensors.

This solution is not very robust to errors in the interference direction.


Subtracting a constant value from all of your samples will just cause all the clipping to happen at -1, and will give you some nasty distortion. The proper way to remove clipping is to normalize the samples.

Keep track of the max value of s1 in your loop, then in a second loop divide all of of your points by the max value, i.e.

a = 0.5;
s1[0] = s[0];
max_s1 = 0;
mean = 0;

for (n = 1; n < N; n++) {
    s1[n] = (s[n-1] * a) + s[n]
    if (abs(s1[n]) > max_s1) {
        max_s1 = abs(s1[n]);
    }
    mean += s1[n];
}
mean /= N;
max_s1 -= mean;
for (n = 0; n < N; n++) {
    s1[n] = (s1[n] - mean) / max_s1;
}

EDIT: For bonus points your can also subtract the mean from your filtered signal prior to doing the normalization. This will ensure you get maximum level without distortion from clipping.


The baud rate/symbol rate is the inverse of the time distance $T_\mathrm{s}$ of two consecutive symbols. It determines the bandwidth of the modulated signal and thus comes into play when simulating frequency dependent components like filters, a wireless channel, a digital-to-analog converter, and so on. If in a simulation system the carrier modulation is simulated in the passband the ratio of baud rate and carrier frequency has to be taken into account. Imagine your carrier is modelled by $\sin(2\pi (f_\mathrm{c}/f_\mathrm{s}) n)$ with a period of $N=f_\mathrm{s}/f_\mathrm{c}$, where $f_\mathrm{c}$ is the carrier frequency and $f_\mathrm{s}$ is the sample rate of the simulation system. In the simulation system, one symbol then has to have the length $N_\mathrm{s}=T_\mathrm{s}f_\mathrm{s}$ or $N_\mathrm{s}=R_\mathrm{s}f_\mathrm{s}$, where $R_\mathrm{s}$ is the baud rate.

If however, you "just have a bit vector converted into a complex vector of symbols" without further upsampling, then the baud rate does not have to be taken into account.

Usually, communication systems are simulated in the equivalent baseband, i.e. the generally relatively fast oscillating carrier is not modelled by a sine wave. Instead, the signal is thought of being shifted to the carrier frequency and only the frequency band of interest that contains the modulated signal is taken into account. Here, the baud rate again determines the width of this band.

If the receiver applies oversampling, i.e. if the sampling rate $f_\mathrm{A}>R_\mathrm{s}$, then the received signal has to be downsampled to reconstruct the transmitted symbols. Oversampling is usually applied for clock recovery or a fractionally spaced equalizer. In a simulation system you could assume perfect sample time and $f_\mathrm{A}=R_\mathrm{s}$. But this depends on what should be simulated, of course.

Update

Generally, in a receiver applying oversampling, the received samples corresponding to a symbol are not identical due to various impairments: noise, inter-symbol-interference, impulse shape, etc. To reconstruct the transmitted symbols, several steps are necessary including clock recovery, synchronization, channel estimation and equalization, and resampling.


Find the Fourier transform of $g(k)$ from $G(z)$ for frequency $=1/2$

First of all it's very important to realize that without any further information it is not possible to uniquely identify the corresponding time-domain sequence from the given $G(z)$. The reason is that you need to know the associated region of convergence (ROC). If the ROC is not known then for the given $G(z)$ there are three possible time-domain sequences. One of them is causal, another one is anti-causal, and the third one is two-sided.

If the ROC is given as $|z|>\frac12$, or, equivalently, if you know that you're looking for a causal sequence, then the inverse transform is uniquely defined. It is similar to the one you found, but there is an important difference:

$$g[k]=\left[\left(-\frac13\right)^k+\left(-\frac12\right)^k\right]u[k]$$

You need the unit step $u[k]$, which makes the sequence causal. The sequence you found extends to infinity for positive as well as for  negative indices $k$, and does not even have a $\mathcal{Z}$-transform.

Now that you know that the ROC includes the unit circle, you know from the answers to this question that in order to find the Fourier transform you can simply evaluate the $\mathcal{Z}$-transform on the unit circle, i.e. set $z=e^{j\omega}$. What they didn't tell you is that this is only possible because the unit circle is included in the ROC. Otherwise you can't do that because the Fourier transform wouldn't even exist.


How to tell whether an image has a low gamma (Finding out imaging system tone map)

You should take a photograph of an object with known color values, such as Macbeth chart:

By using the lowest row of grayscale colors, you will be able to find out the tone-mapping (sometimes incorrectly called gamma) of your imaging system.


Why Fourier Transform doesn't show obvious periodicity?

An FFT doesn't necessarily show a fundamental periodicity, but instead spectral periodicities, which might be higher harmonics (or overtones) of the periodicity that you might perceive visually.

Also, some spectral peaks can be hidden under (the windowing artifacts of) a relatively  large enough DC offset.  If you subtract the mean of the data before doing the FFT, these frequency peaks might show up more clearly.


Estimate a microphone's transfer function

Your $H^{-1}$ is effectively what is called a zero-forcing equalizer for the channel represented by $H$. That channel incorporates DAC, amplifier, speaker, your room, the microphone, amplifier, ADC and several filters necessary to fulfill sampling conditions.

Problem is that for any frequency where $H$ has a notch, $H^{-1}$ has to have a very strong amplification! Sadly, that means that the frequencies where your estimate is based on the least received energy have the largest amplitude in the resulting inverted channel model. That leads to amplified noise.

There's a lot of ways to estimate a channel inverse; the MMSE estimator is pretty popular in communications.

Because you ask for a Python implementation: GNU Radio is usually used from Python, and there's multiple equalizers that come with it (try first!); and, there's modules like gr-adapt which allow for calculating (and using!) channel-reversing filters adaptively. What I especially like about gr-adapt is that its documentation and comments are plentiful and come with references.

For example, have a look at https://github.com/karel/gr-adapt/tree/master/examples :

Here, gr-adapt's LMS Filter is used to estimate (and invert) the channel (which is simulated using the Decimating FIR Filter block).


Basic Confusion About the DFT and Convolution

As I understand your question, it boils down to: why is it ok, to take the convolution of a periodic signal only on an interval rather than over infinity? The DFT or discrete nature of the signal does not really play a role in your question. It merely provides the basis of assumung your two signals $f$ and $g$ as periodic, on the interval $[0,T]$ such that $$f(t+T)=f(t-T)=f(t)$$ and obviously $g(t+T)=g(t-T)=g(t)$.

Now assume that you define the convolution over the whole $(-\infty,\infty)$ like so:
$$y(t):=(f\ast g)(t)=\int_{-\infty}^\infty d\tau\, f(t-\tau)g(\tau)$$

From the periodicity above follows that
$$\int_{-\infty}^\infty d\tau\, f(t-\tau)g(\tau) = \int_{-\infty}^\infty d\tau\, f(t+T-\tau)g(\tau)$$
which you can intepret intuitevely to mean that the convolution should also be periodic. I write "intuitively" because the integral is divergent as we will see below.

Furthermore you can split the integral into a sum of integrals where $\tau$ goes over intervals of length $T$, such as $[n\cdot T,(n+1)\cdot T]$ with $n$ a signed integer number in $(-\infty,\infty)$. Doing that we have (after some simplifications)

$$\begin{eqnarray}
y(t)&=&(f\ast g)(t)=\int_{-\infty}^\infty d\tau\, f(t-\tau)g(\tau) \\
&=&\sum_{n=-\infty}^\infty \int_0^Td\tau\,f(t-\tau)g(\tau)
\end{eqnarray}$$
which would be divergent just because of the sum. Thus, the convolution of periodic functions produces a divergent integral. Sticking to this definition of the convolution would make it impossible to define the convolution like that. At least one of the functions must go to zero sufficiently fast, to produce a finite integral. Your periodic functios, however, will violate this decay condition. Thus it makes sense to window the functions with a rect function that captures all the information in one period.

So you see why it would make sense to define the cyclic convolution as just the integral

$$y(t)=\int_0^Td\tau\,f(t-\tau)g(\tau)$$

and use the periodic extension, since we already provided some intuition as to why setting $y(t+T)=y(t)$ makes sense. Thus the cyclic convolution carries all the information you would have from the regular definition but avoids the divergence due to the infinite sum. I hope this helps.


Enhance image in matlab - exp operation - imageJ software-LPR

You have to learn to debug your programs. Try formula for some data:

>> 255*exp(255/46.01)
  ans =   6.5091e+004

So formula has bug (explanation in link to forum is valid. Read it carefully). Valid formula:

>> exp(255/46.01)
 ans =  255.26

But try calculation with real data from your image:

>> in = (im(10,10)+1)
in = 254
>> exp(in/46.01)
ans =  403.43    % very strange result!!!!
>> in/46.01
ans = 6          %what is it????

Reason is quite simple - type of in is uint8!!!

Valid formula:

result(i,j) = exp(double(in)/46.01);

result:

imshow(uint8(result));


Unstable plant transfer function identification

See the image (found in this article).

Inject a signal into the loop somewhere (I show it below as going between the controller and the plant).  Then measure the input to the plant and the output from the plant.  Note that the image is a bit confusing -- it assumes that the normal command to the loop is 0 or some constant.  It would be better to actually measure right at the output of the plant.

Then use the input to and output from the plant to estimate the plant response.

You can do this for swept sine, impulses, random noise, and steps.  If you're concerned about noise, you can take long samples -- when I do swept-sine measurements using this technique, I set things up so that I can take arbitrarily long measurements at each frequency, to make sure I can capture the real signal in the presence of noise.  If you're doing steps, you can get the response to an arbitrarily high number of steps, average, and then do ARMA or whatever method you wish to extract the plant transfer function.  Ditto random data or impulses.


How to estimate filters using conjugate gradient?

This is a typical deconvolution problem that you can solve either by transforming to the frequency domain where convolution is a simple multiplication:

\begin{equation}
F\left \{ {\bf{I}} \right \}= F\left \{ {\bf{z}} \right \} F\left \{ {\bf{f}} \right \},
\end{equation}

where ${\mathcal{F}}$ denotes the Fourier transform (or DFT) so $f$ will be:

\begin{equation}
{\bf{f}}={\mathcal{F}}^{-1}\left \{ \frac{{\mathcal{F}}\left \{ {\bf{I}} \right \}}{{\mathcal{F}}\left \{ {\bf{z}} \right \}} \right \}
\end{equation},
where ${\mathcal{F}}^{-1}$ is the inverse Fourier transform.

Otherwise, you can also represent in a matrix form, using a Circulant (kind of Toeplitz)  for the convolution operation. The frequency domain calculation though does not involve inverting (or pseudo inverting) a matrix.

If there's additive noise in your problem, please see Wiener deconvolution (both time and frequency domain implementation)


Under certain conditions, a $T$-periodic function can be represented by its Fourier series

$$x(t)=\sum_{k=-\infty}^{\infty}a_ke^{jk\omega_0t}\tag{1}$$

with $\omega_0=2\pi/T$. The function $x(t)$ can be complex-valued, i.e. have non-zero real and imaginary parts. Note that generally the Fourier coefficients $a_k$ are also complex-valued (even for real-valued $x(t)$).

Now if $x(t)$ is real-valued, i.e., $x(t)=x^*(t)$, you get

$$x(t)=x^*(t)=\sum_{k=-\infty}^{\infty}a^*_ke^{-jk\omega_0t}=
\sum_{k=-\infty}^{\infty}a^*_{-k}e^{jk\omega_0t}\tag{2}$$

Comparing $(2)$ with $(1)$ you see that for real-valued $x(t)$ you get the condition $a_k=a^*_{-k}$, i.e. the Fourier coefficients show Hermitian symmetry. This condition is not satisfied for general complex-valued functions $x(t)$. However, the representation $(1)$ is valid no matter if $x(t)$ is real-valued or complex-valued.


Matlab video processing of Contracting tissues. how to analyze movement?

You could try a simple non-linear denoising filter (3x3 or 5x5 median filter) on your image sequence to make your detection more robust to noise ; but I get the feeling that the moving details on these images or quite thin and that too much filtering would harm your detection. And it looks like you are already thresholding the images before comparison, so the only noise that matters here are pixels moving above and below the threshold at each frame.

To get a beating rate you could also compute the auto-correlation of the sequence. It will show a secondary maximum at the beating rate. The autocorrelation is actually an interesting representation for your task because:

A very regular beat will mean a very sharp secondary peak on the AC sequence ; while a more irregular beat will show a flattened, spread out secondary peak in the AC sequence.
The ratio of the autocorrelation secondary peak to the first peak ($r_0$) can be used as a good measure of beat strength.

But back in time domain you could also consider:

The kurtosis of the sequence as a measure of "peakedness" (and thus beat strength). Other measures of peakedness include the ratio of geometric mean to arithmetic mean.
The standard deviation of the time intervals between peaks as a measure of beat regularity. The downside is that you will have to detect peaks, which will be a bit imprecise in your case especially given the low temporal resolution of your data (and this is why I thought about doing things from the AC function first!)


What is the unit of autocorrelation function?

It depends.  Normalized autocorrelations have no units because the units are divided out as part of the normalization process.

Non-normalized autocorrelations have the original data's units squared.