Kalman filter : simple code example I read lots of things about Kalman filtering, but in order to fully understand it, I would probably need to see it working on some data.

Would you have a minimal example (Python code or any other language) showing what it does on some real data $x[n]$, where $n$ is the time?

I created a simple example with Scilab/XCOS. It is a simulation of lowpass filter with noise. It is observed by a kalman filter. You can also insert some uncertainties in the system model. With XCOS you can simulate the system. The results will be plotted automatically.
Here is the Downloadlink.

Little help with scilab:
Go to the directory with standard unix command

cd /home/workingdir/Kalman_Example

There you should execute getd() to load all functions (.sci-files) in the directory.
Then you can execute the modeling.sce  file, which is just an initializing script:

exec('modelling.sce')

Then typexcos simulation_kalm1.zcos into the console. A graphic window will open. Press the play button and the simulation will start.


Given some parameters, how do you determine how much of an image sensor an object will consume? If I have the details of a scene and a camera, how can I tell how much of my camera's image sensor will be taken by a certain object in space?

For example, if I have a painting that is 1m by 1m, it is straight ahead 10 m away from my camera.  If my image sensor is 5 mm by 5 mm, and I fix the focal length of the camera, how do I determine how much of the image sensor will be taken up by the painting?

The basic idea is that the ratio of the distance/focal length will be the same as the physical/image ratio.

So, assuming a focal length of 100mm, this gives:

$mag=10m/100mm=100$
$size_{image}=1m/mag=10mm$

Note, the distance is the distance measured from the focal point, which will be 1 focal length from the sensor, so that might factor into the equation as well.

See this calculator to play around with things a bit more.


Counting blinks of a LED I'm working on an algorithm that counts blinks of a LED. I have series of non-neagtive integers. They represent brightness of a LED in time. I give you two examples of these integer series below. Visually it's very obvious that LED blinked 7 times in example 1 and 8 times in example 2. But how can this be calculated by a computer. Any ideas? I don't even know if I'm asking this question on the right site. I'm new to this stuff but any help would be appreaciated.

Example 1:
Data of Example 1

Example 2:
Data of Example 2

You can do something like this in Python (setting data to your original signals):

import numpy as np
data = data1

max_val = np.max(data)
thresh = max_val/2
thresh_crossings = np.where(np.diff(np.signbit(np.array(data) - max_val/2)))[0]
print(len(thresh_crossings)/2)

This code returns 7 for your first signal and 8 for the second one.

First I find the maximum value. Then I arbitrarily set a threshold (in this case, half that value), taking advantage of the fact that the signals are pretty clean and there are not any spurious peaks. Doing so leads to something like this:

Then you just have to find how many times the signal crosses that threshold, and divide it by two (because it will cross it twice for each peak, once when it rises and afterwards when it falls).

If you have any additional information a priori of the frequency of the peaks, you could add some hysteresis to make the algorithm more robust (i.e. if there are small zigzags in the signal, avoid counting them as threshold crossings).


What is the difference between a sampled time system and a discrete time system? I am trying to evaluate whether I should control a system in simulink using sampled data system or discrete time system.

The discrete model is as follows:

The sampled time system is as follows

What is the functional difference between these two and which approach is preferred in realistic modelling and control?

The first one is a purely discrete model. You need to select a discrete solver in the simulation options to accurately simulate this system. The second system however is actually continuous because of the holders, so you need a continuous time solver.

The preferred approach depends on what your system is. If you have a discrete model to work on, you should prefer the first one. If you have an actual physical system that you want to control with a microcontroller, you should prefer the second one.


Maths behind why removing values from DFT is not a brick wall filter I sort of understand that when removing values from a signals DFT is not the same as a brick wall filter since it will only zero the frequencies that are perfectly periodic. Other spectrum frequencies who are not exactly periodic will have their energies spread over the whole DFT result. Is this right?

I'd like to know how you can mathematically prove that zeroing frequencies is not a brick wall filter using the DFT and inverse DFT.

A true brickwall filter would require an infinite impulse response in the time domain, since the inverse Fourier Transform of a rectangular window (brickwall frequency response) is a Sinc function. Since for the DFT the number of samples in the frequency domain is equal to the number of samples in the time domain, this would require you to have an infinite number of samples in the frequency domain and then require you to null an infinite number of samples.

For anything less, the impulse response in time becomes an aliased Sinc function which no longer corresponds to a true brickwall frequency response.

If you restrict your frequency space to just the DFT samples however, you will of course achieve a brickwall response following the rules that all samples outside of your rectangular window (in frequency) are zero. However in the implementation of a filter we are often concerned with the frequency response of signals that exist on a continuous frequency domain, in which case as described above you will not be able to achieve a perfect brickwall frequency response.


Why is the total signal response response of zero input + zero state, $y(t) = y_0(t) + h(t)\star x(t)$, not an LTI? For the zero input + zero state response in continuous time linearly time-invariant systems, why is the $y(t)$ equation not "technically" considered an LTI? I read this in a journal and there was no justification. How would you change the equation to make it a true LTI?

The response of the system is given by

$$y(t)=h(t)\star x(t)+y_0(t)\tag{1}$$

where $x(t)$ is the input signal, $h(t)$ is some well-behaved function, $y_0(t)$ is a function which is independent of $x(t)$, and $\star$ denotes convolution. The model $(1)$ can be used for describing systems with non-zero initial conditions that would be LTI were it not for the non-zero initial conditions (such as an initially charged capacitor).

For $y_0(t)=0$, $h(t)$ would be the impulse response of the corresponding LTI system. For $y_0(t)\neq 0$, $h(t)$ is not the impulse response, neither can the system be characterized by an impulse response.

For $y_0(t)\neq 0$, the system described by $(1)$ cannot be linear, which is easily shown by considering a response $y_1(t)$ to an input signal $x_1(t)$. For a linear system, the response to a scaled version of $x_1(t)$ must be a scaled version of $y_1(t)$, with the same scaling factor. So the response to $ax_1(t)$ must be $ay_1(t)$ if the system is linear. But from $(1)$ we have

$$y_1(t)=h(t)\star x_1(t)+y_0(t)$$

and for the scaled input we get

$$y(t)=h(t)\star ax_1(t)+y_0(t)=ay_1(t)+(1-a)y_0(t)\neq ay_1(t)$$

for $y_0(t)\neq 0$ and $a\neq 1$. Consequently, the system described by $(1)$ is not linear.

In a similar way it can be shown that the system is also not time-invariant. The reason is again the function $y_0(t)$, which is independent of the input signal $x(t)$. If $x(t)$ is shifted, the output doesn't shift accordingly, because $y_0(t)$ doesn't shift.


Calculate sum of absolute difference per pixel I am trying to implement a research paper just to verify my results. I have shared the code (written by me) and the description of the process given in the research paper. Unfortunately my results are a bit different from his results that means I am doing something wrong. Kindly guide me in the right direction.

Here is the link to the paper: A parametric model for perceptual video quality estimation.

I have tried for mobile video sequence which can be downloaded from the following link: Yuv video sequences.

He gets 3.600 as SAD/pixel value while I am getting around 1.

Code:

fileName = 'mobile_cif.yuv';
    width=352;
    height=288;
    numFrames=30;
    SAD=0;

    %for loop to traverse & process from frame '1' to 'last' frames
    for t = 1 : numFrames
    if t < numFrames %(Need to traverse to n-1 frames)
    result=loadFileYuv(fileName,width,height,t);
    currFrame = result.cdata;    %reading Current frame
    result=loadFileYuv(fileName,width,height,t+1);
    nextFrame = result.cdata;    %reading Next frame
    % First, take the absolute value of the difference at each pixel
    myAbsDiff = abs(double(currFrame) - double(nextFrame));
    % Then, sum over all pixels
    out = sum(myAbsDiff(:));
    out=out/(width*height); %(Taking Average with respect to pixel)
    out=out/64;%(Dividing By Block size to normalize)
    disp(out);
    SAD=SAD+out;
    else
        continue;
    end

    end
    disp(SAD/numFrames);  %(Taking Average with respect to Frames)

Process Description

Sum of absolute difference (SAD) is a simple video quality metric used for block comparison and for moving vectors calculations. Each frame is divided into small blocks (i.e. $8 \times 8$ pixels) and for every block in one frame the most similar (minimum SAD) block in next frame is find. This minimum sum of absolutes differences is assigned as the SAD for each block in each frame (up to the $n − 1$ frame). Then all the SAD values are averaged for each frame and for all the frames in the clip, and divided by the block area, for normalization.

I think there are multiple things to consider here:

1) Where are you dividing the image into 8x8 blocks? If you don't divide like this, then normalization by 64 is meaningless, and you are just computing SAD over all pixels. Maybe I'm seeing something wrong? You should do the normalized SAD in those small windows and then sum and normalize. And since 3600 is not a normalized value, I don't think that he normalizes.

2) Have you considered that your function for loading the video or the codec you use is actually giving you slightly different values than the one he uses? If your input is different, you will obviously get different results.


Find "base" frequency of pulse wave I have some samples that are read from a square pulse wave and I'm looking for a way to find the characteristics of this wave.

It's a pulse wave, let $T_1$ be the non-pulse duration and $T_2$ to be the pulse duration, usually $T_1$ is bigger than $T_2$.

I would like to know the base frequency, e.g., the frequency that includes $T_1+T_2$.
Can I do this using an FFT?

Of course. Sample a few cycles, window, and do a FFT. The highest peak will be found at the base frequency (T1+T2) of your wave. With "spectral reassignment" you might get more precise results.


Filter design with a phase delay specification? Given a desired filter magnitude response, plus acceptable linear phase and minimum phase FIR filter approximations, how can one design a filter with the given response and a delay specification of some value between minimum phase and linear phase?

Least square errors works well with FIR filters. General IIR filters are more difficult and typically require an iterative search algorithm. One specific type of IIR, named warped FIR filters, can also match arbitrary amplitude and phase response with a least square errors approach.


Phase response of $H(f)=e^{-j2{\pi}ft_0}$ Given is the impulse response: $$h(t)=\delta(t-t_0)$$.

I calculated $$H(f)=e^{-j2 \pi ft_0}=|H(f)|\cdot e^{j\varphi(f)}$$.

Now, the magnitude response of $H(f)$ is:

$$\begin{align}
|H(f)| &=\sqrt{\Re\{H(f)\}^2+\Im\{H(f)\}^2} \\
 &= \sqrt{\cos^2(2{\pi}ft_0)+\sin^2(2{\pi}ft_0)} \\
 &= 1 \\
\end{align}$$.

Using: $$e^{-j \theta}=\cos(\theta)- j \cdot \sin(\theta)$$ and $$\cos^2(\theta)+\sin^2(\theta) = 1$$.

But now I wonder about the solution of the phase response. I know that the solution is:

$$\varphi(f)=-2 \pi f t_0$$.

I would be thankful if someone could show me how the phase response was calculated, since following

$$\varphi(f) = \arctan \left(\frac{\Im\{H(f)\}}{\Re\{H(f)\}} \right) = \arctan \left(\frac{-\sin(2 \pi f t_0)}{\cos(2 \pi f t_0)} \right)$$

I instead obtain $$\varphi(f)=-\arctan\big(\tan(2{\pi}ft_0)\big)$$ which differs from the given sample solution.

Wouldn't it be much easier to just compare the left-hand side and the right-hand side of

$$e^{-j2\pi ft_0}=|H(f)|e^{j\phi(f)}\tag{1}$$

to see that $|H(f)|=1$ and $\phi(f)=-2\pi ft_0$?


How to determine the central frequency from a vectorial spectrum analyzer? I have taken measures of an antenna with a vectorial spectrum analyzer. I need to obtain the central frequency at which the antenna works. In order to do so, I have selected S11 (reflection) and I have plotted the "log magnitude" of the spectrum, I obtained something like this:

I know that the central frequency is obtained from one of the peaks shown on the plot. But in this case, I obtain three peaks. Which one should I take if they are pretty much the same? (no peak is considerably greater than the others)

The S11 represents the relative amount of signal that is reflected from the antenna, so the lower the better, because lower means more radiated by the antenna.

Usually, we consider the antenna to be "working" at a given frequency if its S11 is lower than -10 dB (of course, this requirement can be changed according to the application. spatial antennas are more around -20 dB for example).

Therefore, your antenna seems to be working at three distinct frequencies (with the s11 < -10 dB criterion). Is this a problem ? It depends on the hardware you put behind it during operation.


Computational Complexity of Polyphase Resampling does someone know the approximate computational complexity of a polyphase resampling algorithm?

I am trying to find a good argument to prefer a polyphase resampling of Scipy over the FFT resampling method. In its documentation for resample_poly() it says:

  This polyphase method will likely be faster than the Fourier method in scipy.signal.resample when the number of samples is large and prime, or when the number of samples is large and up and down share a large greatest common denominator.

My data meets those criteria. I resample 15 channels * 30k-50k samples from 10kHz to 300Hz and practice shows it is way faster than the FFT method.

Yet I would like to have a referable source on those criteria.

I found a few papers on the computational complexity of FFT algorithms to be around O(n*log(n)) with n=n_channel*n_samples.

Rational Resampling

10 kHz -> 300 Hz is a rational resampling with relatively benign factors:

$$ \frac{300\,\text{Hz}}{10\,\text{kHz}}=\frac{3\cdot10^2}{10^4}=\frac3{100}\text,$$

meaning that you'd interpolate by 3 and decimate by 100. That's exactly what you're achieving using resample_poly:

Polyphase Decimators

Remember that for decimation by 100 without aliasing, you need to filter your signal to $\frac1{100}$ of it's original bandwidth, which might take quite a lot of filter taps. Sadly, you throw away 99 of 100 samples that you filter in the following decimation step.

Using polyphase techniques, you can cleverly rearrange the $\frac1{100}$-band filter taps so that you don't even do these 99 calculations – and reduce the effort of your decimation filter by the decimation factor, here, 100.

This basically happens by taking your long filter and "deinterleaving" it to 100 disjunct sets, forming 100 polyphase components of your original filter. If you put in the input samples in all of these partial filters, and sum up the output samples of these component filters correctly, you get exactly the same filter. It's really just a way of "arranging" your filter differently. However, if you know which output samples you're going to throw away, you can omit pushing in each input sample to 99 of 100 polyphase filter components – so that you instead "round-robin" the input samples to each of the polyphase filter components. That saves 99% of calculations!

Polyphase resamplers

But what about the upsampling by 3?
Well, another beautiful aspect of polyphase rational resamplers is that you can show that the order of "interpolate by 3 first, eliminate images by $\frac13$-band filter, then apply $\frac1{100}$-band filter to avoid aliasing, then throw away 99 of 100 samples" can be reordered to decimate-filter-interpolate. That's awesome because:

the long filter runs at the lowest rate in the system
you can omit the shorter filter, because it's far more relaxed than the longer one, anyway.

Why is that faster than the $\mathcal O(N+Nlog(N))$ of an FFT filter?

Aside from the aforementioned dramatic reduction in number of necessary operations, remember that even for the shorter polyphase component filters, you can apply FFT filters, if they are long enough so that this pays; in my experience, for rates as benign as yours, that's usually not the case. Compare a post by Tom Rondeau about comparison of FFT FIR filters to polyphase FIRs.

Generally, the beauty of FFT filters lies in the fact that you can reduce the number of operations you need to do to convolve a piece of signal with a filter tap vector. On the other hand, you buy that algorithmic elegance with totally non-linear memory accesses when computing that FFT. Now, if your filter is long, so must be your FFT, and your CPU can't keep all of both the signal vector and the tap vector in localized caches. Rule of thumb says that a single fetch from RAM is worth 200 to 500 Floating point operations[citation needed] on modern x86 computers.

However, even if the amount of operations you need to handle is quite a lot larger in the non-FFT FIR filter, your memory access is (if your FIR implementation isn't totally braindead) inherently linear – and your CPU and its memory controller are really good at fetching linear memory, and even pre-fetching it. So, that works out in your favor quite a long time, until the complexity advantage of the FFT filter really beats the crap out of the naive FIR.

Another aspect of that linearity is easier vectorization using SIMD instructions – if you want to feel disappointed by scipy, install GNU Radio (easy if you're on current Linux distros or OS X, still doable on windows), run volk_profile (takes a while to find the optimum SIMD implementations for your machine), and then use the following flow graph to measure resampling throughput:

It does 5 million samples per second for me. Scipy's not getting quite close to that.


Interpreting dataset of a power spectrum I am currently working with kaldi, and are trying to make spectrogram plots of different audio files. Kaldi is a speech recognition toolkit.

Kaldi frames its input in frames of 25 ms with a hop_length of 10 ms.

Kaldi have some code that makes spectograms of an audio files, but not much information is given on how things are stored...
The sample rate used for the audio files is 16 KHz.

An example of dataset could be fash-b-an251:

which is an audio file of duration 1 sec.
the matrix has 98 rows, and 257 columns.

Plotted:

A different dataset could be fash-b-an253

which is an audio file of duration duration: 0.70s
The matrix has 69 rows and 257 columns.
Plotted:

257 must be the size of the fft length, since it is consistent.
What is the rows - how it is related to the length of the audio file.
And can the length of audio file be determined using number of rows(Probably), and if,then how?..

I am also a bit unsure about the unit of the values in the plots..
Is is a log power spectrum, or just a power spectrum?

Spectogram of the First audio file with matlab:

So the unit are not the same, what is it stored as??

1 second at 10ms "hops" gives 100 hops, which correspond to the rows.

In DSP parlance, the frame (the 25ms of samples) is called window size, and is the sample interval that is fed into the DFT (FFT). These samples are then multiplied by a window function prior to calculating the DFT.

25ms with $f_s = 16 \rm{kHz}$ gives 400 samples per frame. From your data, the DFT is performed with a vector size of 512, so probably the frame is then 0-padded from 400 to 512 samples, then fed to the DFT.

The DFT will then produce its output as a vector of frequency samples of size 512, where each frequency corresponds to multiples of $f_s/512 = 31.25 \rm{Hz}$. So your first column corresponds to frequency 0 Hz, then second column to 31.25 Hz, third to 62.5 Hz, etc.

Only half of the output is shown; since your signal is real, the other half of the DFT output (from 8000 to 16000 Hz) will be symmetrical to the first half. Also, since the sampling frequency is 16 kHz, then you cannot represent signals beyond 8 kHz.

Why 257 columns? Probably because the last column corresponds to 8 kHz (inclusive), so you have 256 frecuency intervals, therefore 257 frequency samples.

The plot looks like log of the magnitude of the DFT. The exact units you can deduce studying your audio signal. Try yourself with Octave / Matlab, using the spectrogram functionality there.

The DFT of 1 frame gives you 1 row of your spectrogram. Step forward 10ms, get another frame, repeat.

P.S. Sometimes spectrograms are easier to visualize with time on the horizontal axis.


How Do I Measure the Time Duration of a Finite Length Discrete Sequence? Assume I have a five-sample time-domain sequence (none of the five samples are zero valued) and the time period between each pair of samples is one second. Measured in seconds, what is the time duration of my five-sample sequence?

The following graphic should clear up why a useful answer is 5 seconds, consistent with conversions between the discrete-time and continuous-time domains. A good example of this are discrete time approximations to continuous time integration or differentiation where the inclusion of $\Delta T$ (sample period) is made. With all such approximations, the time duration attributed to each sample is one sample interval. For the Forward Euler approximation specifically, each sample covers the time duration up to and not including the next sample. Such as in the first graphic the time duration for a sample at $t=0$ is one sample interval of one second; up to and not including the next sample at $t=1$ when sampled at one sample per second, with a total time duration for the one sample as one second.
The bottom figure then shows how this would apply to five samples starting at $t=0$ resulting in a total time duration of 5 seconds.

The plot below shows the Forward Euler discrete-time approximation for the integration of a continuous time waveform $x(t)$, which would be given as follows in the formula below and demonstrated in the graphic, best showing how each sample carries the time duration of the sample period (and as David has pointed out in the comments, this equally represents the zero-order hold approach in the output of a typical Digital to Analog converter):
$$\int_0^{NT}x(t)dt \approx \sum_{n=0}^{N-1}x[n]T$$

The actual range of time used for each sample is specific to the analog-digital translation approach used, but would cover one sample duration in all cases, leading to a total duration of 5 seconds for the OP's case.


About the orthogonality of columns in Zadoff-Chu matrix As we know that Zadoff-Chu matrix is similar to Walsh-Hadamard matrix where every columns in those matrices is orthogonal with the any other column.

For Walsh-Hadmard matrix, the orthogonality is clear and can be demonstrated easily. But regarding the Zadoff-Chu matrix, I don't see that orthogonality is right as shown below:

a = (sqrt(2) + 1j*sqrt(2))/2;
W = [a 1 a -1; -1 a 1 a ; a -1 a 1; 1 a -1 a];  %%Zadoff-Chu matrix

C1 = W(:,1);     %% Take the first column

orth_results = [];
for j = 1 : size(W,2)
    results = (abs(sum(W(:,j).*C1)).^2);  %Project the first column to every column in matrix W
    orth_results = [orth_results results];
end

Following the condition of orthogonality between columns, orth_results should be equals to 0 except in the fist value, but what I get is orth_results = [8.0000         0    8.0000      0];

What's the problem of that? Is there any issue in the matrix itself ?

For complex matrices, the concept of orthogonality is replaced by unitarity. If $^H$ denotes the  conjugate transpose of a matrix, we should check that $W$ is a unitary matrix (up to a scaling factor) via equations:
$$W^*  W = WW^*  =\lambda I$$
with $\lambda\neq 0$, which you can test positively with:
W'*W
W*W'

What is missing in your loop, and especially in the complex scalar product, is the conjugation on one of the complex vectors. This is explained for instance in Wikipedia/Dot product for Complex vectors.
You could for instance replace C1 by conj(C1).


Creating highpass from lowpass in Python I want to create high pass filter from low pass filter in Python. My code:

h_lowpass = lp_design_window(fc, N, window)
dirac_delta = np.concatenate((np.ones(1), np.zeros(h_lowpass.shape[0] - 1)))
h_highpass = dirac_delta - h_lowpass

where lp_design_window is the function with arguments: cutoff frequency, number of coefficients and window function - and returns coefficients of low pass filter using window method. It attenuates desired frequencies, but not enough. Is it something wrong with dirac_delta? Or maybe I should search for the problem somewhere else?

You need to line up your direc_delta with the center tap of your low pass filter. Put the 1 in the middle at N/2.


How to circularly shift a signal by a fraction of a sample? The shift theorem says:

  Multiplying $x_n$ by a linear phase $e^{\frac{2\pi i}{N}n m}$ for some integer m corresponds to a circular shift of the output $X_k$: $X_k$ is replaced by $X_{k-m}$, where the subscript is interpreted modulo N (i.e., periodically).

Ok, that works fine:

plot a

N = 9
k = [0, 1, 2, 3, 4, 5, 6, 7, 8]
plot ifft(fft(a)*exp(-1j*2*pi*3*k/N))

It shifted by 3 samples, as I expected.

I thought you could also do this to shift by fractions of a sample, but when I try it, my signal becomes imaginary and not at all like the original:

plot real(ifft(fft(a)*exp(-1j*2*pi*3.5*k/N)))
plot imag(ifft(fft(a)*exp(-1j*2*pi*3.5*k/N))), 'b--'

I didn't expect this at all.  Isn't this equivalent to convolving with a real impulse that's been shifted by 3.5 samples?  So the impulse should still be real, and the result should still be real?  And it should have more or less the same shape as the original, but sinc interpolated?

If you want the shifted output of the IFFT to be real, the phase twist/rotation in the frequency domain has to be conjugate symmetric, as well as the data.  This can be accomplished by adding an appropriate offset to your complex exp()'s exponent, for the given phase slope, so that the phase of the upper (or negative) half, modulo 2 Pi, mirrors the lower half in the FFT aperture.  The complex exponential shift function can also be made conjugate symmetric by indexing it from -N/2 to N/2 with a phase of zero at index 0.

It just so happens that the appropriate offset for phase twists or spirals, that complete an exact integer multiples of 2 Pi rotations in aperture, to be conjugate symmetric in aperture, is zero.

With a conjugate symmetric phase twist vector, the result should then end up as a circular Sinc interpolation for non-integer shifts.

Elaboration by OP:

Your choice of k = [0, 1, 2, 3, 4, 5, 6, 7, 8] is producing an asymmetrical complex exponential:

If you use k = [0, 1, 2, 3, 4, -4, -3, -2, -1] instead, you get a Hermite-symmetric complex exponential:

plot(fftshift(exp(-1j * 2*pi * 0.5/N * k)))

and now when you use the same exponential formula to shift by 0.5 or 3.5 samples, you get a real result:

plot ifft(fft(a)*exp(-1j * 2 * pi * 0.5/N *k))
plot ifft(fft(a)*exp(-1j * 2 * pi * 3.5/N *k))


FM modulation audio 50 Hz – 16 kHz MATLAB I'm trying to modulate in FM an audio signal $x(t)$ coming from a microphone which works from 50 Hz to 16 kHz in MATLAB.

But the problem is that I want a FM Bandwidth of 25 kHz at the output of modulator.

I have some problems to do this:

The FM bandwidth for any signal using Carson's rule is:

$$BT = \left(\frac{f_d \max\vert x(t)\vert}{W}+ 1\right)W$$

where $W$ is the maximum bandwidth (15.950 kHz in this case) of the modulating signal and $\max\vert x(t) \vert $ is the max value of the modulating signal in $Vp$.

The thing is I have to adjust the deviation sensibility $f_d$ (slope) to get the desired FM bandwidth, but how to know $\max\vert x(t)\vert$ in the microphone datasheet? or even using an amplifier between microphone and modulator?

From a practical point of view, the maximum value that your microphone will send you is limited by the input ADC clipping level. The clipping level is probably stated in terms of the analog reference voltages that the audio codec uses. Moreover, the analog value that's converted into digital may be further normalized or scaled by the audio codec firmware before sent to the calling application.

In the old PC WAV audio standard, 8-bit channels sent a maximum value of integer 255 with a bias of 128. The 16-bit channels used a signed format with maximum corresponding to $2^{15}$. Note that this value is typically mapped to some range and converted to floating point values for further processing such as done in MATLAB. So be aware of that.

The standard DSP diagnosis procedure recommends a test with a known or predictable input signal to observe the mathematical range of the resulting input samples.


Optimal trade-off between oversampling and filter length For some sampling-frequency-preserving operations on Nyquist–Shannon sampled signals, such as:

a shift a.k.a. translation, and
differentiation by applying a derivative filter a.k.a. gradient filter,

the ideal filter for each task has an infinite number of taps with weights obtained by sampling the sinc function or its derivative. In practice, finite-impulse-response (FIR) approximations are typically used.

For a given maximum acceptable level of approximation error, the required minimum number of taps in the approximating FIR filter is reduced if there is an overhead $f_B\ldots f_s/2$ between the highest frequency $f_B$ of the signal and half the sampling frequency $f_s$. From filter design perspective, this overhead obtained by oversampling the input signal by a factor $\beta = f_s/(2f_B) > 1$ creates a transition band of non-zero width and relaxes requirements on the filter, as exemplified by the magnitude frequency response of a simple approximative derivative filter decipted in Fig. 1.

Figure 1. Arbitrarily normalized magnitude frequency response of an ideal derivative filter (diagonal red line) and its approximation (black curve) with impulse response $-3/16,$ $31/32,$ $0,$ $-31/32,$ $3/16$. With ideal oversampling, the frequency response above the highest frequency of the signal (vertical blue line) does not contribute to the approximation error.

However, assuming that the output signal is sampled at the same sampling frequency $f_s$ as the input signal, increasing $f_s$ increases the number of output samples, which increases the number of times a dot product between the $N$ filter weights and $N$ signal samples needs to be calculated. In general, for a 1-d signal, $\beta N$ times as many multiply–accumulate operations (MACs) are needed compared to the number of samples in critically sampled (not oversampled) representation of the input and output signals.

For such a task, which oversampling factor $\beta$ minimizes $\beta N$, requiring the least total MACs in FIR filtering without exceeding a given maximum acceptable approximation error measured by a reasonable error metric?

Upsampling to obtain the oversampled input signal shall be considered ideal and not contributing to the MAC count, as if the input signal was obtained perfectly pre-oversampled.

1-d MMSE derivative filter for uniform spectral prior

With some simplifying assumptions about the signal distribution, a filter with least square frequency response error is the Bayesian minimum mean square error (MMSE) filter for a uniform prior of the frequency spectrum. For such 1-d filters, an oversampling factor $\beta \approx 2\ldots2.5$ seems optimal, as it is in the sweet spot of the cost-benefit curves (Fig. 1).

Figure 1. Mean square frequency response error of least square optimal derivative filters as function of computational complexity as measured by $\beta N$ which is proportional to the total number of multiply-accumulates (MACs), for various oversampling factors $\beta$ for which the filters have been optimized. Also shown are traces for different numbers $N$ of non-zero filter coefficients. No advantage was taken from impulse response anti-symmetry or simple coefficients such as $1$.

More data is needed before general conclusions can be drawn.


Amplitude calculation of the original signal from its FFT I have a signal whose FFT has a one prominent harmonic ($f_1$) and few other harmonics with a lower amplitude. In order to calculate the amplitude of the original signal I can use a well known theorem: a value of $f_1$ divided by half number of data points in the original signal.

But what if I have harmonics $f_2, f_3, \ldots$? Should I take their amplitudes into account for such calculation?

Whoa, clarification is needed. The peak amplitude ($A_p$) of an $N$-length $f_1$ input sinusoid can only be found using

$A_p = 2M/N$

(where $M$ is the the magnitude of the $f_1$ spectral component) if $N$ samples of that $f_1$ input sinusoid represent an exact integer number of cycles.

If that condition holds then you can measure $A_p$ for the $f_1$ sinusoid as well as the peak amplitudes of the harmonics. But those peak amplitude values apply only to the individual sinusoids that comprise your FFT input signal. You CANNOT use the individual measured peak amplitudes to measure the peak amplitude of your nonsinusoidal FFT input.  That was the point made in jojek's comment.


PSD (Power spectral density) explanation I'm trying to understand how the PSD is calculated.  I've looked in a few of my Communication Engineering textbooks but to no avail.  I've also looked online.  Wikipedia seems to have the best explanation; however, I get lost at the part where they decide to make the CDF (Cumulative Distrubution Function) and then for some reason decide to related that to the autocorrelation function.

I guess what I don't understand is, how does autocorrelation having anything to do with calculating the PSD?  I would've thought that the PSD simple be the Fourier Transform of $P(t)$ (where $P(t)$ is the power of the signal with respect to time).

You are right, PSD has to do with calculating the Fourier Transform of the power of the signal and guess what.....it does. But first let's look at the mathematical relationship between the PSD and the autocorrelation function.

Notations:

Fourier Transform:
$$ \mathcal{F}[  x(t)] = X(\omega) = \int_{-\infty}^{\infty} x(t)e^{-j\omega t}dt $$
(Time) Auto-Correlation Function:
$$ R(\tau) = x(\tau) * x(-\tau) = \int_{-\infty}^{\infty} x(t)x(t + \tau)dt $$

Let's prove that the Fourier Transform of the Auto-Correlation function does indeed equal to the Power Spectral density of our stochastic signal signal $x(t)$.

$$ \mathcal{F}[  R(\tau)]= \int_{-\infty}^{\infty} R(\tau)e^{-j\omega \tau}d\tau $$
$$   = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} x(t)x(t + \tau) e^{-j\omega \tau} dtd\tau $$
$$   = \int_{-\infty}^{\infty}x(t) \underbrace{\int_{-\infty}^{\infty} x(t + \tau) e^{-j\omega \tau} d\tau}_{\mathcal{F}[x(t + \tau) ] =  X(\omega)e^{j\omega t}} dt $$
$$   = X(\omega) {\int_{-\infty}^{\infty}x(t)e^{j\omega t}  dt}    $$

$$   = X(\omega)X^*(\omega)= |X(\omega)|^2 $$

What does it all mean?
Note: This explanation is a bit "hacky". But here it goes

The Fourier transform tells us the spectral components of a signal. In our case, the signal is Stochastic; So, trying to calculate the spectral components of the signal will be pointless because, for every realisation of the random process, you will have different expressions for $ \mathcal{F}[x(t)]$.

What if you take the Expected Value of the Fourier transform then? This wouldn't work. Let's take a zero mean signal for example.

$$ \mathbb{E}\{ \mathcal{F}[x(t)] \} = \mathcal{F}[\mathbb{E}\{ x(t) \}] = 0$$

Instead, what if you take the Fourier transform of the square of the signal.
$$ \mathbb{E}\{ \mathcal{F}[x^2(t)] \} = \mathcal{F}[\underbrace{\mathbb{E}\{ x^2(t) \}}_{\text{Av. Power of the Signal}}] $$

The autocorrelation function is essentially the $P(t) $ which you were alluding to.

References:

[1] Communications 1, P-L. Dragotti, Imperial College London

[2] White Noise and Estimation, F. Tobar [Unpublished Report]


Difference in frequency spectrum of continuous time signal and discrete time signal? I am studying chapter 4 of book signal processing first

There i came across a concept which says that(what i have understood)the spectrum of continuous time signal has limited number of lines while spectrum of discrete time signal has infinite number of lines

Why?what is the reason for this difference?

https://imgur.com/a/6IENWha

above link
Is snap of  frequency spectrum diagrams of both continous time and discrete time signals

below link is relevant theory of book
https://imgur.com/a/mH7L7D6

This is the classical sampling theorem. When you sample a continuous signal $x(t)$, you are basically multiplying it by a sample train $s(t) = \sum_{-\infty}^{+\infty}\delta(t-kT_s)$, the value at consecutive $T_s$ being your samples. In frequency domain, it effect is to convolve $X(\omega)$ and $S(\omega)$

$$S(\omega) = \sum_{-\infty}^{+\infty}\delta(\omega-2\pi k/T_s)$$

So the resulting spectrum would be
$$
X_s(\omega) = \sum_{-\infty}^{+\infty}X(\omega-2\pi k/T_s)
$$
There are infinite shifted copies of $X(\omega)$. To begin with your original signal was a single-frequency signal which traced the circumference with frequency $f_m$. When you sampled it at $f_s$ rate, you got infinite copies at $f_m \pm kf_s$, where $k$ is any integer in $(-\infty,+\infty)$ hence it referred to as infinite number of lines after sampling (discrete time).


Feature Selection by Filtering so there is this paper I'm reading and trying to understand fully: Towards Practical Identification of HF RFID Devices
http://dl.acm.org/citation.cfm?id=2240276.2240278

I don't want to link the PDF here, but maybe you can catch it from ACM or somewhere else.

Now what we need to know: They record a signal for $5 \mu s$ with 4 GS/s. From the normalized amplitude of the signal $l$ at time $t$ called $f(t,l)$ they obtain $F(\omega, l)$ by FT.

$$
F(\omega, l) = \frac{1}{\sqrt{M}} \sum_{m=0}^{M-1} f(t,l) \exp (-2 \pi i \frac{t\omega}{M})
$$
Where M is the length of signal considered and $0 \leq t \leq M − 1$ is time. Now quoting:

  We then remove from $F(\omega, l)$ its DC component and the redundant part of the spectrum;
  we denote the remaining part of the spectrum by $\vec{s}_{l}$.

So now for the spectrum between 0 to 100 MHz the vector $\vec{s}_{l}$ contains amplitudes for 20000 frequencies $(D=20000)$. To reduce the dimensions the authors filter the signal with a bandpass.

  A projected vector $\vec{b}_{l}$, also called filtered features, is extracted
  from the basic features $\vec{s}_{l}$ using a bandpass filter transformation $W_{BPF}$:
  $$
\vec{b}_{l} = W_{BPF}^{t} \vec{s}_{l}
$$
  The feature extraction from N captured responses for a given RFID device is
  then given by $G = W_{BPF}^{t}S$ where $G$ is an array of  $\vec{b}_{l}$ and
  $S = [\vec{s}_{0} \dots \vec{s}_{l} \dots \vec{s}_{N}]$ is a matrix that
  contains the basic features of N captured responses.

  In our implementation, the transformation $W_{BPF}$ consisted of a Chebyshev Type I bandpass filter with filter order of 100 [Williams and Taylors 1988]. The filter was used to keep the frequencies in a target frequency range enclosed between the upper and lower bandpass edge frequencies, denoted by Fp1 and Fp2 respectively; the latter were experimentally determined. We chose the filter order to be high for better bandpass filtering quality. Using this technique the basic feature dimensionality can be effectively reduced to only tens of frequencies (dimensions) comprised between Fp1 and Fp2.

Later on in the paper they state the following:

  Therefore, we applied feature selection by filtering with parameters Fp1 = 11 MHz and Fp2 = 15 MHz. This procedure reduced the features (fingerprint) dimensionality to D = 40.

And now my question (lately).

Why is the dimensionality now 40? If we have amplitudes of 20000 frequencies between 0-100 MHz and filter from 11 to 15 MHz shouldn't there be 800 left?
How does this bandpass matrix $W_{BPF}^{t}$ look like?

Thanks for reading and greetings,

para

The problem is that the following statement- "So now for the spectrum between 0 to 100 MHz the vector s⃗ l contains amplitudes for 20000 frequencies (D=20000)"- is faulty.

A signal that lasts $5\mu s$ is captured at 4 GS/s.  That is a total of 5e-6 * 4e9 = 20000 samples.  The Nyquist frequency is half the sample rate, or 2 GHz.  They take the Fourier transform of the samples ($F(w,l)$).  At this point, assuming the samples are real (as opposed to complex), there are 20000 frequency bins representing -2 GHz to +2 GHz.  The negative frequency information is redundant since it mirrors the positive frequencies.  I believe that when they said "We then remove from F(ω,l) its DC component and the redundant part of the spectrum" they were talking about getting rid of the 0 Hz and negative frequency information.

That would leave you with 9999 frequency bins representing 200 kHz to 2 GHz.  To make life easier, though, let's bring back the DC bin and say that we have 10000 frequency bins representing 0 Hz to 2 GHz.

Hopefully it is clear now why "So now for the spectrum between 0 to 100 MHz the vector s⃗ l contains amplitudes for 20000 frequencies (D=20000)" is faulty.  There are 10000 bins and they represent 0 to 2 GHz, not 0 to 100 MHz.  They now bandpass the signal down to 4 MHz bandwidth.  $(4e6 / 2e9) * 10000 = 20$ frequency bins.  I believe they get 40 dimensions from those 20 bins by counting both the amplitude and phase of the bins as features, $20 * 2 = 40$.


Fourier Transform of shifted Impulse I'm new here and I hope you can help me.
I started studying the Fourier Transform now at University and I have a lot of doubts about this subject.
My teacher sent to us one list of exercises and I don't know how I can solve the Fourier Transform of this function:

$$
x(t) = \delta \left(\frac{t-T}{a}\right)
$$

I don't know how to use the properties of scaling and shifting in this case, it would be great if you could explain to me the steps to solve this question.
Thanks for your help :)

Technically, the impulse $\delta(t)$ is called a distribution, and not a function, but for the purposes of your first course in Fourier transforms,  what you need to know is that $\delta(t)$ has the sifting property

  $$\displaystyle\int_{-\infty}^\infty x(t)\delta(t) \,\mathrm dt = x(0) ~\text{provided that }x(t) ~\text{is continuous at }t=0 \tag 1$$

which sifts the single value $x(0)$ out the function $x(t)$ and discards all other values. Also, any integral involving impulses can be manipulated via the change of variable technique into an integral of the form above to get the value of the integral. For example,

\begin{align}
\int_{-\infty}^\infty x(t)\delta(t-T) \,\mathrm dt
&= \int_{-\infty}^\infty x(\tau+T)\delta(\tau) \,\mathrm d\tau
&\scriptstyle{\text{set }\tau = t-T, ~ dt = d\tau}\\
&= x(T)
\end{align}

and so we can think of $\delta(t-T)$ as a time-shifted impulse that plucks the value $x(T)$ (instead of $x(0)$) out of the function $x(t)$.
Similarly, for $a >0$,
\begin{align}
\int_{-\infty}^\infty x(t)\delta(at) \,\mathrm dt
&= \int_{-\infty}^\infty x\left(\frac{\tau}{a}\right)\delta(\tau) \,\frac 1a \mathrm d\tau
&\scriptstyle{\text{set }\tau = at, ~ dt = \frac 1a d\tau}\\
&= \frac 1a x(0)
\end{align}
Figure out why we can think of $\delta(at)$ as the magnitude-scaled impulse
$\frac 1{|a|}\delta(t)$ and then apply both of the above ideas to
evaluating the Fourier transform integral
$$\int_{-\infty}^\infty e^{-j\omega t}\delta\left(\frac{t-T}{a}\right)
\, \mathrm dt$$


Dynamic Time Warping - Comparing Values Ok, so I'm trying to compare two different speech signals and I have come into a problem. Here goes:

I have split the signal into blocks, and I have computed the MFCC coefficients of each block. I then use a DTW algorithm to compare the (inputted) signal to the training signal.

Here is the data:

The training data is the MFCC values of a signal of someone saying "No". If I have the inputted signal of someone saying no the values then give:

  13.9462,
  55.8784,
  38.3383,
  29.9468,
  32.7136,
  24.8893,
  34.0734,
  24.3645,
  39.329,
  20.4847,
  31.1939,
  29.8841,
  25.7655,
  28.222,
  23.1643,
  33.4366,
  .....,
  ......

Where the max data for this dataset is: 97.4834

Here is the data for someone saying a different word:

  28.6696,
  65.8777,
  44.2725,
  31.6083,
  42.6541,
  38.4104,
  26.6311,
  34.9188,
  37.2065,
  25.2479,
  41.5969,
  54.2681,
  37.0685,
  26.2073,
  33.9836,
  38.7847,
  28.3622,
  67.8788,
  74.9075,
  .....,
  .....
  Where the max value in this dataset is: 97.5609

My guess (reading through the algorithm) is that the smallest difference should be the correct result. However, I do not know if it is possible to either:

1) Calculate the MAX value for each dataset and then see which value is smaller

2) Get the first variable and compare it.

Does anyone have any suggestions to how I would compare these kind of values?

Thank you

The result of DTW is a matrix $N\times M$, where $N$ is the size of the first MFCC-dataset in comparison and $M$ is the size of the second one. Now to make a decision of how close two datasets are to each other you need to see at the normalized bottom-right value value of the DTW matrix. If it's less than some pre-defined threshold, then you can say that two datasets are of a similar nature. In order to compute the threshold you need to run DTW multiple times and plot a ROC curve. The following sample DTW implementation in Matlab should help:

function cost = dtw (A, B)
  [rowsA colsA] = size (A);
  [rowsB colsB] = size (B);

  D = inf (rowsA + 1, rowsB + 1);
  D(1) = 0;

  for i = 2:rowsA + 1
    for j = 2:rowsB + 1
      d = norm ( A(i - 1, :) - B(j - 1, :) );
      D(i, j) = d + min ([ D(i - 1, j - 1) D(i - 1, j) D(i, j - 1) ]);
    endfor
  endfor

  cost = D(end) / (rowsA + rowsB);
endfunction


Programming Wavelets for Audio Identification How exactly is a wavelet used digitally?

Wikipedia states:

  For example, a wavelet could be created to have a frequency of Middle
  C and a short duration of roughly a 32nd note. If this wavelet was to
  be convolved with a signal created from the recording of a song, then
  the resulting signal would be useful for determining when the Middle C
  note was being played in the song.

Would this be a data structure holding e.g {sampleNumber, frequency} pairs?

If a wavelet is an array of these pairs, how is it applied to the audio data?

How does this wavelet apply to the analysis when using an FFT?

What is actually being compared to identify the signal?

Convolution can be used to find the note. Its best shown visually .

  Would this be a data structure holding e.g {sampleNumber, frequency} pairs?

Yes, the signal would be. The wavelet can be stored as a function but it can also be {sampleNumber, frequency} pairs.

  If a wavelet is an array of these pairs, how is it applied to the audio data?

Convolution:

$f*g = \int f(x)g(t-x)dx$

Convolution is commutative so it does not matter whether f or g is wavelet or signal or vice versa.

  How does this wavelet apply to the analysis when using an FFT?

In some applications in signal processing, convolution is calculated using the convolution theorem.

$\mathcal{F}(f*g) = \mathcal{F}(f) \mathcal{F}(g)$

so the convolution of f and g can alternatively be calculated by

$f*g = \mathcal F ^{-1}[ \mathcal{F}(f) \mathcal{F}(g)]$

  What is actually being compared to identify the signal?

The signal is being compared to the wavelet

  How exactly is a wavelet used digitally?

There are many applications, besides note identification, they can be used in lossy signal compression. See the Dirac codec.

Cross Correlation is a similar operation which can be used for identifying wavelets


phase response of an analog filter I am trying to understand the relation between phase response of an analog filter and its influence on the quality of a signal with PSK modulation.

For example I want to design a BPF for HF communication system.

The pass band is [1.5 .. 30]MHz. The Phase response is non-linear as far as I understand. Will it damage the phase modulated signal, say PSK-8 ?

Which of two - phase response or group delay is more important to be analysed?

Any references where I can read about that?

-Thanks in advance

A non-linear phase response is equivalent to group delay variation, and similar to multi-path fading can cause inter-symbol interference if it is severe enough, but either can be corrected with channel equalization. In this case it would be a static correction so would not be overly complicated to implement in the transmitter, receiver or both. Factoring the equalization into the cascade of two equal filters and then implementing it in both the transmitter and receiver would also provide the benefit of SNR optimization (against white noise) as a matched filter (but the actual performance benefit of doing that may be negligible compared to the added complication, so just throwing out the thought).

The metric to analyze would be group delay variation since that would show it directly. A linear phase component is only a delay and does not provide any distortion, so the group delay which is the derivative of the phase with respect to frequency would convert any linear phase to a constant value and highlight the variation which is the destructive component. One approach is to simulate the filter response and analyze the eye diagram to assess the signal quality and if further equalization and or changes to the analog filter design would be needed.

The graphic below depicts this for the case of multipath, but multipath is simply a filter with phase (and amplitude) variation across the signal bandwidth.


Why do we need to write the noise as randn(2,1)./sqrt(2*SNR) instead of randn(1,1)/sqrt(2)? I saw this code, about the error probability of on/off signal ,on the internet

Ts=10^4;
SNR_dB=[-3 0 3 6 9];
Ne=zeros(1,length(SNR_dB));

for i=1:length(SNR_dB)
    SNR=10^(SNR_dB(i)/10);
    for t=1:Ts
        s=((sign(rand()-0.5))+1)/2; %create the "on"signal or "off" signal ."on"signal=1,"off" signal=0
        r=randn(2,1)./sqrt(2.*SNR); %it seems that it is noise
        if (s+r)>0.5
            s_hat=1;
        else
            s_hat=0;
        end
        if abs(s_hat-s)>0
            Ne(i)=Ne(i)+1;
        end
        y=Ne./Ts;
    end
end
Pe=zeros(1,5);
for i=1:5
    SNR=10^(SNR_dB(i)/10);
    Pe(i)=erfc(sqrt(SNR/2)/(2^0.5))/2;
end
semilogy(SNR_dB,y,'-r',SNR_dB,Pe,'-k')

I don't understand why can the noise be written as

  r=randn(2,1)./sqrt(2.*SNR)

,why do we need to write as randn(2,1)./sqrt(2*SNR) instead of randn(1,1)/sqrt(2)?Is anyone know about this?

Not sure if this will help, but let me show how to generate a noisy signal of a desired SNR, assuming additive white Gaussian noise (AWGN) model with fixed signal power.

First the basic definition of SNR (signal to noise ratio) is :
$$ \text{SNR} = \frac{\sigma_x^2}{\sigma_v^2} $$

where ${\sigma_x^2}$ is the signal power based on a WSS random signal model, and ${\sigma_v^2}$ is therefore the noise power (which I take here as white Gaussian noise).

SNR can also be defined in decibells:
$$ \text{SNR}_{dB} = 10 \log_{10} \left( \frac{\sigma_x^2}{\sigma_v^2} \right) = 10 \log_{10}( \text{SNR}) $$

Assuming your signal power $\sigma_x^2$ is fixed, then you can adjust the noise power $\sigma_v^2$ to achieve a desired SNR with your AWGN signal model. The necessary noise power is found from required SNR as:

$$ \sigma_v^2 = \frac{ \sigma_x^2 }{ SNR } $$

Now, Matlab's built in Gaussian function randn() generates standard normal with zero mean and unit variance. Then to get a desired variance (hence noise power) you should scale it with the square root of the target power as shown:

sigx2  = 5.1;          % signal variance (power)
SNR    = 23 ;          % SNR = sigx^2 / sigv^2
sigv2  = sigx2 / SNR;  % noise variance (power)
N      = 1024;         % number of noise samples
v = sqrt( sigv2 ) * randn(N,1); % White Gaussian noise

The reason why you multiply with the square root of the power (or variance) is based on the following property of the variance operator:

$$ \text{Var}(X) = \sigma_x^2  \implies \text{Var}(aX) = a^2 \sigma_x^2$$


Z-transform of not quite an upsampler I know the z-transform of an upsampler is:

$$
y[n] = \begin{cases}
x(n/L) &n=0,\pm L, \pm 2L, ...\\
0&otherwise
\end{cases}
\longrightarrow  Y(z)= X(z^{L})
$$

if $x[n]_L$ is defined to zero for $n <0$ and $n \ge L$ and to be non-zero over interval $[0, L-1]$...then, what's the Z-transform of:

$$
y[n] = \begin{cases}
x[n]_L & 0 \le n \le L-1\\
0&otherwise
\end{cases}
\longrightarrow  Y(z) = ???
$$

its kind of like:

$$
y[n] = x[n]\ u[-n + L]\ u[n]
$$

I have this z-transform for LPF that is similar, but doesn't match:

$$
\begin{cases}
y[n]=a^{n}&0 \le n \le L-1 \\
\\
0& otherwise
\end{cases}
\longrightarrow  Y(z) = \frac{1-a^Lz^{-L}}{1-a{\ z}^{-1}}
$$

Ok, its a trick question...its $Y(z) = X(z)$ because $y[n]$ and $x[n]_L$ are identical sequences.

$$
y[n] = \begin{cases}
x[n]_L & 0 \le n \le L-1\\
0&otherwise
\end{cases}
\longrightarrow  Y(z) = X(z)
$$


How does the assumption that symbols are equi-probable hold If I draw a number uniformly between zero and one, what is the probability that they are equal? Mathematically, it should be zero but I don't recall why? Can somebody please help in explaining why the probability should be zero?

For example, there is still a chance,  that when generating six random numbers they will be all 1 or all 0? So if the probability is zero then why in digital communications do we assume that the symbols are $iid$ and that if the symbols are 0 and 1, 50% of them will be 1 and the rest will be 0 i.e., why do we assume equi-probable symbol probabilities when theoretically they cannot be equal?
In Matlab, I tried my best to generate a large sequence of 0/1 but never did I get equi-probable symbol probabilities.

Am I misunderstanding some thing here?

If I draw a number uniformly between zero and one, what is the
  probability that they are equal? Mathematically, it should be zero but
  I don't recall why? Can somebody please help in explaining why the
  probability should be zero?

Avoiding formal definitions of (Lebesgue) probability measure, an informal way is thinking the probability at a point of a uniformly distributed random variable is $\frac{1}{\textrm{number of points}}$. Then the number of points of real values between zero and one is infinity, thus the probability must be $\frac{1}{\infty} = 0$.

  For example, there is still a chance,  that when generating six random
  numbers they will be all 1 or all 0? So if the probability is zero
  then why in digital communications do we assume that the symbols are
  $iid$ and that if the symbols are 0 and 1, 50% of them will be 1 and
  the rest will be 0 i.e., why do we assume equi-probable symbol
  probabilities when theoretically they cannot be equal?  In Matlab, I
  tried my best to generate a large sequence of 0/1 but never did I get
  equi-probable symbol probabilities.

The previous arguments change if we divide the interval $[0,1]$ to several sub-interval of width $\Delta$. Thus the probability falling on each sub-interval is $\frac{1}{1/\Delta}=\Delta \neq 0$.

In your example, $\Delta = 0.5$ yields two sub-intervals and each has probability $0.5$. That's where the equiprobable symbol assumption holds.

Experimentally, e.g. with Matlab, you (almost surely) never get exactly the probability $\Delta$ because the number of tests is always finite.


Ideal low pass filter output at given sampling frequency Consider the signal $\cos(30t)$ sampled at $w_s=40 rad/s$ using a unit impulse train. The sampled signal is filtered with an ideal low pass filter with unity gain and cutoff frequency $w_c = 40rad/s$. Find the resulting output.

The solution is given as $\frac{20}{\pi}(cos(30t)+cos(10t))$

Attempt:

Taking the Fourier transform of $cos(30t)$ gives $\frac{\pi}{2} (\delta(\omega-30)+\delta(\omega+30))$.

Using the unit impulse train:
$X(\omega)_{T_s} = \frac{1}{T_s}\sum X(\omega-k\omega_s)$ but this gives $\frac{1}{T_s}\sum \frac{\pi}{2} (\delta(\omega-70)+\delta(\omega+10))$.

If this is an acceptable approach, it's unclear where the cosine arguments and factor of $\frac{20}{\pi}$ is coming from. Help is appreciated.

That $\frac{1}{T_s}$ factor in your sampled signal expression is $\frac{20}{\pi}$.

$$\omega_s = 40 \ rad/sec$$
$$2\pi f_s = 40$$
$$f_s = \frac{40}{2\pi} = \frac{20}{\pi}$$

The sampled signal is basically given by :
$$X_{sampled}(f) = f_s \sum^{\infty}_{k=-\infty} X(f - kf_s)$$
Meaning : Sampling a signal with impulse train at sampling frequency $f_s$ gives you, images of $X(f)$ scaled by $f_s$ and centered at $kf_s$.

You can check this answer to understand why : Sampling with Impulse Train explained pictorially

The original signal is at frequency $\frac{15}{\pi}$ and being sampled at sampling frequency $\frac{20}{\pi}$. And hence the frequency components in the sampled signal will be : $\pm k\frac{20}{\pi} \pm \frac{15}{\pi}$.

Notice that for $k=0$, you get the original signal and for $k=\pm 1$, you have an alias term $\cos{(2\pi (\frac{20}{\pi} - \frac{15}{\pi})t)} = cos(10t)$ which will get inside the lowpass filter cutoff frequency.

And, that will give you the final result as :
$$f_s (\cos{(2\pi \frac{15}{\pi}t)} + \cos{(2\pi \frac{5}{\pi}t)})$$
$$ = \frac{20}{\pi} (\cos(30t) + \cos(10t))$$


Comparison between 2 images in Frequency domain As you can see there are 2 images from the Tom and Jerry cartoon program.

In the 1st image both Tom and Jerry are present. But in the 2nd one,only Tom is present.Now,we can clearly see this difference in spatial domain. But i want to compare these images in Frequency domain.

So can anybody explain how do we compare these 2 images and find out the missing part from the 2nd image i.e. image of Jerry using Fourier transform ?
Also,is it possible to find out location  of Jerry in the image using any frequency transform method? If yes can you explain how?

1st question: there are some simple steps to change the spatial domain to freq domain, exp adding some zero to avoid ISI, you can read the procedure from a source book, you can use fft2 command in matlab

2nd question: Surely you can not see the differences in freq domain as clear as spatial domain, because to determine the value of one pixel in freq domain all the pixels in the spatial domain will be effective


Feature matching of images without corners I am trying to stitch aerial images using feature-based registration. Most image pairs are matched well, such as this:

Others are not matched due to lack of common features:

As you can see, the images contain almost no corners (only fuzzy edges).

Is there a way to register such images? Please note that I need homography so template matching is not possible. The images are taken by flying drone wiggling in the air, hence the images have perspective distortions and feature matching is therefore essential.

Please note I have tried many point features (Harris, SIFT, SURF, MOPS) with similar (bad) results.

My ideas:

Use some line-based (Hough transform?) or area-based features (MSER?)
Use more images or video and then recover homography from dense data using Horn-Schunck or Lucas-Kanade

Here you can download the source images - it contains two pairs of images, the easier pair and hard-to-align pair:

Download original images

I have referred to Stage I.D of this tutorial. Hope this helps.
http://www.robots.ox.ac.uk/~vgg/practicals/instance-recognition/index.html#stage-id-improving-sift-matching-using-a-geometric-transformation

When the features have scale and orientation assigned (e.g. SIFT features have these properties), you can compute similarity transform between each corresponding feature pair. The similarity transform has 4 parameters: Tx, Ty - translation, R - rotation and S - scale.

The tentative correspondences are then filtered using RANSAC-like algorithm such that only correspondences consistent with common similarity transform (up to a predefined threshold) are used.


Frequency of the wave in frequency domain If we have a 1-dimensional wave in time domain, it can be represented in frequency domain with x axis indicating the frequency of the wave and y axis indicating amplitude/magnitude of the wave. But what about a 2-dimensional wave in case of images in Frequency Domain? Because in 2-dimensional wave we have a cordinate (x,y) to denote the magnitude. So what is the frequency in this case?

If you're working with a 2 dimensional data set, you need to do a 2D FFT on it (which can be dome using many 1D FFTs: doing all the rows and then all the columns, or if you prefer, the columns and then the rows).

The 2D FFT gives you a 2D array of complex numbers.  Horizontal waves will give peaks on the horizontal axis.  Vertical ones will give peaks on the vertical axis, and diagonal waves will give peaks away from both axes.


Matlab: White noise with flat/constant power spectrum I am relatively new to signal processing techniques and Matlab and need a bunch of test data in the form of white noise, as defined on the Wikipedia page, with constant, flat power spectrum.

I tried the matlab built-ins wgn(m,n,p) and randn(m,n) to get vectors of white noise, and then used different combinations of the periodogram/pwelch/dspdata.psd functions to check what I believe to be the power spectrum. I cannot say I am 100% certain on the basis of documentation and other forum sourcing, so do correct me if I am mistaken.

The resulting plots I get do not look even close to a flat spectrum but instead inherit the "zick zacky" curve of the original random data. The question I inevitably had to ask myself: Do you even know what to expect and how to understand the results of the functions and how a flat power spectrum looks like? I don't, hence I am here, looking for clarification.

I am looking for answers which go beyond the Wikipedia pages etc., as they left me confused.

You need to average a bunch of FFTs of white noise to approach the average power spectral density.

see this and this.


Taking inverse Fourier transform in column-wise and solve it in row-wise $\DeclareMathOperator{\FFT}{FFT}\DeclareMathOperator{\IFFT}{IFFT}$Assuming I have a matrix $X$ of size $64\times16$. Taking the $\IFFT$ for it in column-wise, I means that $Y = \IFFT(X)$;

Is it possible to get a relationship between every row in $X$ and its corresponding row in $Y$?

No. Your problem doesn't put any constraints on the columns, so they're totally independent. Then, transforming them to a different base doesn't change that, at all.

If the entries in different columns had nothing to do with other columns, that won't change. The fact that you arranged your data in a matrix doesn't magically allow any operation to introduce relationships that weren't there.


What are integration scale and differentiation scale? In scale-adapted Harris detector, the scale adapted second moment matrix is defined by:
$$\mu(x, \sigma_I, \sigma_D) = \sigma_D^2\ g(\sigma_I) *\left[ \begin{array}{cc}
L_x^2(x, \sigma_D) &L_xL_y(x, \sigma_D) \\
L_xL_y(x, \sigma_D) &L_y^2(x, \sigma_D)
\end{array} \right]$$,
while $\sigma_I$ is the integration scale, and $\sigma_D$ is the differentiation scale. What are they? Why use them?

UPDATE

According to @Dima's answer, I think $\mu(x,\sigma_I, \sigma_D)$ is actually something like this,
$$\mu(x,\sigma_I,\sigma_D) = \sigma_D^2\cdot g(\sigma_I)*\left(g(\sigma_D)*\begin{bmatrix}I_x^2(x) &I_xI_y(x) \\
I_xI_y(x) &I_y^2(x)\end{bmatrix}\right)$$, where $I_x,I_y$ are the 1st order derivatives. Is it right? If so, why bother convolving twice? And why multiply $\sigma_D^2$?

The $\sigma_I$ determines the scale level at which the Harris corners are computed. Coarser scales (higher values of $\sigma_I$) correspond to larger corners.

The $\sigma_D$ is effectively the window size, over which the derivatives are summed to generate the entries of the matrix. If $\sigma_D$ is too small, then the detector will be seriously affected by noise.  If it is too large, then small corners may wash out.  In the original non-scale-adapted Harris detector you only have $\sigma_D$.

By the way, this is the first time I see the terms "integration scale" and "differentiation scale", and I am not sure I like them. The "differentiation scale" is actually the Gaussian envelope, over which the derivatives are integrated. :)


Determining the 'amount' of speaking in a video I'm working on a project to transcribe lecture videos. We are currently just using humans to do the transcriptions as we believe it is easier to transcribe than editing ASR, especially for technical subjects (not the point of my question, though I'd love any input on this). From our experiences we've found that after about 10 minutes of transcribing we get anxious or lose focus. Thus we have been splitting videos into ~5-7 minute chunks based on logical breaks in the lecture content. However, we've found that the start of a lecture (at least for the class we are piloting) often has more talking than later on, which often has time where the students are talking among themselves about a question. I was thinking that we could do signal processing to determine the rough amount of speaking throughout the video. The idea is to break the video into segments based on the amount of typing they require, as opposed to splitting based on time.

I've done a little research into this, but everything seems to be a bit overkill for what I'm trying to do. The videos for this course, though we'd like to generalize, contain basically just the lecturer with some occasional feedback and distant student voices. So can I just simply look at the waveform and roughly use the spots containing audio over some threshold to determine when the lecturer is speaking? Or is an ML approach really necessary to quantify the lecturer's speaking?

Hope that made sense, I can clarify anything if necessary.

Appreciate the help as I have no experience with signal processing.

The most straightforward and easy thing you could do is rectify the signal (i.e. take the absolute value of its samples) and then run a moving average filter to get an approximate estimate of the "envelope" of your signal. Imagine the moving average filter like estimating the mean of the signal values between samples 0-100, then 1-101 then 2-102 and so on.

The result of this will be an inaudible slow varying waveform. You can then do a histogram on its values and use that to estimate a threshold which could then be used to "tag" the interesting parts of your video.

However, I feel that this is a lot of hard work for this kind of problem.

Therefore, I would also like to suggest that you use a video editor that presents the video "signal" as two separate ribbons, one for the actual video and another representing the amplitude of the sound. You can then quickly inspect the audio ribbon and estimate quickly the parts that could be seen as "interesting" and use the video editor to either directly transcribe those areas or cut the video down into smaller segments. Unfortunately, at the moment, I can`t recall a freely available video editor with this feature but if I do I will amend this response.


Should you scale the FFT bins by $1/N$ where $N$ is the number of points in a transient signal? I am going over the following tutorial. In section 2.1 the author says

  To calculate the $N$ point FFT the Matlab algorithm 1 can be used. Here, after taking the FFT, its magnitude is calculated and the bins are scaled by $1/N$.

Now the author doesn't explain as to why the bins are scaled by $1/N$. At first glance it looks like some kind of normalisation -- but why?

I did some further research and it seems there different ways to define the power in a signal, depending on what you are looking at. See slide 3 of this.

  http://www.hep.ucl.ac.uk/~rjn/saltStuff/fftNormalisation.pdf

Yes, some kind of normalization needs to be done and it's a matter of convention which one to use. To explain, let $\mathbf{F}$ be a DFT matrix containing ${\rm e}^{-\jmath 2\pi \frac{mn}{N}}$. Then you can show that $\mathbf{F} \cdot \mathbf{F}^H = N \cdot \mathbf{I}_N$. This factor $N$ needs to be accounted for. Now, you can do either of the following:

Define DFT via $\mathbf{D} = \mathbf{F} \cdot \mathbf{d} \quad \rightarrow\quad$ IDFT becomes $\mathbf{d} = \frac 1N \mathbf{F}^H \cdot \mathbf{D}$.
Define DFT via $\mathbf{D} = \frac 1N \mathbf{F} \cdot \mathbf{d} \quad \rightarrow\quad$ IDFT becomes $\mathbf{d} =  \mathbf{F}^H \cdot \mathbf{D}$.
Define DFT via $\mathbf{D} = \frac{1}{\sqrt{N}} \mathbf{F} \cdot \mathbf{d} \quad \rightarrow\quad$ IDFT becomes $\mathbf{d} = \frac{1}{\sqrt{N}} \mathbf{F}^H \cdot \mathbf{D}$.

All three of them can be used and none of them is more or less correct than the other. Matlab uses the second one, Maple the third (if I remember correctly).

If you need the units to resemble Fourier integrals, normalizing by $\Delta t$ or $\Delta f$ can make sense as well, as you found in the slide set. Depends on what you want to do with the result really.


Introduction to IP : Laplacian of Gaussian? I am very new to the concept of computer vision (and image processing), and trying to understand the algorithm used for edge detection.
One thing I'm currently struggling to figure out is the Laplacian of Gaussian

In this case, how to you determine the value of x and y? I thought the centre of the kernel was the origin $((x,y) = (0,0))$, but it doesn't seems like the right number.

Also, if there are any resources that you would recommend to go through, I would really appreciate them.

You are right about the $x$ and $y$ values. The center of the matrix is $(0,0)$ and the corner points are $(\pm 4,\pm 4)$. But they obviously wanted integer values in the matrix, so they simply scaled the LoG function by a factor of $482.75$ (just to get a decent range). Evaluating the function with this scale factor gives you (for the lower right quarter of the total matrix):

LoG =

  -40.000000  -23.086993    0.294243    5.218347    2.080704
  -23.086993  -11.762411    3.077873    4.839310    1.745668
    0.294243    3.077873    5.409023    3.361998    0.998325
    5.218347    4.839310    3.361998    1.456020    0.365518
    2.080704    1.745668    0.998325    0.365518    0.081641

Rounding should give you the final result. However, if you look closely, you'll see that it doesn't (at least not for all $(x,y)$ values). If you check the matrix on this page, you'll see that for the same $\sigma$ it is also different. So people make mistakes and they just copy from each other.


Intuitive explanation of subspace There are many techniques in signal processing that use eigen analysis (MUSIC, SVD, eigen decomposition, etc) that result in signal and noise subspaces.The mathematical definitions for signal subspaces are abundant, but what is the intuitive, tangible explanation of what a subspace is representing? More importantly, how does one interpret the values of a subspace? What exactly does this result provide and what is an example of how one would use it? Nearly any topic I can think of in signal processing has very intuitive explanations of complex topics - but I've yet to see a good one for subspaces.
EDIT: The crux of the question is, what is the intuitive explanation of subspace as it applies to signal processing algorithms and applications (i.e., not the linear algebra explanation)?

TL;DR:  Subspaces are low-dimensional, linear portions of the entire signal space that are expected to contain (or be close to) a large part of the observable and useful signals  or transformations thereof, with additional tools that allow us to compute interesting things on the data
We are  given a set of data. To manipulate them more easily, it is common to embed them, or represented them in a well-adapted mathematical structure (from the plenty of structures  we have in algebra or geometry), to perform operations, prove things, develop algorithms, etc. For instance in channel coding, group or ring structures can be better adapted. In a domain called mathematical morphology, one uses lattices.
Here, for standard signals or images, we often suppose a linear structure: signals can be weighted, added: $\alpha x+ \beta y$. This is the base for linear systems, like traditional windowing, filtering (convolution), differentiating, etc.
So, a mathematical  structure of choice lies in vector spaces. Vector spaces equipped with tools: a dot product (that can be used to compare data), a norm (to mesure distances). These tools help us compute. Indeed, energy minimization and linearity are strongly related.
Then, a data of $N$ samples naturally lives in the classical linear space of $N$ dimension. It is quite big (think of million-pixel images). It contains an awful lot of other "uninteresting" data: any $N$ dimensional "random" vector. Most of them are and will  never be observed, have no meaning, etc.
The reasonable quantity of  signals that you can record, up to variations, is very small relatively to the big space.
Even more, we are often interested in structured information. So if you subtract noise effects, unimportant variations, the proportion of useful signals is very very tiny within the whole potential signal space.
One very useful hypothesis (heuristic, to help discover) is that those interesting signals live close together, or at least along regions of the space that "make sense". An example: suppose that some extraterrestrial intelligence has no other detection system than  a very precise dog detector. They will get, across the Solar system, almost nothing, except many points located on something vaguely looking like a sphere, with large empty spaces (oceans), and  sometimes very concentrated (urban areas). And the point cloud moves around a center, with a constant periodicity, and rotating on itself. Those aliens have discovered something!
Anyway, the partial-sphere looking point cloud is interpretable... maybe a planet?
So, our dog point cloud could have been fully 3D, but they are concentrated on a 2D surface (lower dimension), that seems relatively regular (in altitude) and smooth: most dogs live at intermediate altitudes.
These smooth low-dimensional parts of space are sometimes called smooth manifolds or varieties. Their structure and operators allow to compute things. For instance: distances, distributions, etc. Inter-dog distances make more sense when computed along the Earth surface (in spherical 2D coordinates) than directly through the planet with the standard 3D norm! But this  can still be complicated to deal with. Let us simplify this a bit more.
Looking a little closer, the dog points are almost located on close-to-flat surfaces: countries, even continents. Those flat surfaces are portions of linear (or affine) subspaces. Still, you can now compute inter-dog distance, more easily, and design an algorithm for dog matching that will make you rich.
The story continues a bit. Sometimes, natural data does not assemble around a clear structure, directly. Unveiling this inherent structure is at the core of DSP. To help us in this direction, we can resort to data transformations to concentrate them better (Fourier, time-frequency, wavelets), filtering.
And if we find a suitable subspace, most algorithms become simpler, more tractable, and so on: adaptive filtering, denoising, matching.
[ADDITION] A typical use is the following: a signal can be better concentrated with a well-chosen orthogonal transform. In the meantime, a zero-mean random Gaussian noise remains Gaussian under an orthogonal transformation. Typically, the covariance matrix can be diagonalized. If you sort the eigenvalues in decreasing ordre, the smallest ones tend to flatten (they correspond to noise), and the highest more or less correspond to the signal. Hence, by thresholding the eigenvalues, it because possible to remove the noise.


Blind Source Separation of Real World Mixture I have been toying with both non-negative matrix factorization (NMF), independent component analysis (ICA) and independent vector analysis (IVA) for separating speech mixtures. I'm trying to separate two sources using recordings from two omni-directional microphones. When I generate a mixture synthetically  in an instantaneous fashion all algorithms provide excellent separation. When I generate convolutive mixtures using room impulse response data from the web the performance is still good. However all the algorithms fail miserably when the input mixture is a real-world recording, like a male and a female reading aloud in a room.

All the examples I could find online, for example this one, also use synthetic convolution. This paper seems to suggest the same, when it says

  In controlled laboratory experiments, BSSD algorithms perform very well. In the real world,
  however, they flounder. If we focus our research more on the practical acoustic engineering necessary
  to obtain signals with less reverberation and noise, then our current algorithms will thrive.

I'd like to know if any new algorithms have been discovered that deal with this kind real-world mixture and provide satisfactory results.

Independent Vector Analysis which extends the idea of ICA seems to be the most robust blind source separation algorithm. It's able to separate real world convolutive recordings into surprisingly high quality components. IVA deals with the permutation problem by modeling each source as a multivariate random variable. For more details the original IVA paper by Taesu Kim can be found here: Independent Vector Analysis: An Extension of ICA to Multivariate Components.


Is this image good enough for OCR? Hello I'm trying to ocr on a special image.

I've done adaptive threshold , polartocartesian transformation.

this is my original image

this is adaptive threshold

this is polartocartesian transformation

what do you think is this image good enough for OCR or do I have to optimize my preprocessing steps ?

It looks like using the polartocartesian transformation is destroying the characters (look at the "M", "E" before and after...). I think it will be a better idea to think about some rotations mechanism using imrotate.

Any way, I'd try the following additional pre-proc:
1. you might want to perform some morphological operations on the image after the adaptive threshold.
2. I would also try labeling the image (bwlabel) and removing too large objects (most of the background noise in your image is connected so it will be a huge object - in terms of pixel counting)

Have fun!


Using k Nearest Neighbour & Dynamic Time Warping for keyword spotting As stated in the title, I am currently researching a method using k-NN, specifically, as either an alternative or as a supplement to the DTW algorithm in keyword spotting based on MFCC. I have read through various answers on this forum (mainly Speech recognition using MFCC and DTW(Dynamic Time Warping)?) and I would like to know if anyone could direct me to somewhere I can find more research or implementations which employ the k-NN.

If you want to use the k-NN in conjunction with DTW, then it's quite easy. Assuming that you have the speech corpora of keywords you should do the following at the training stage:

Each keyword should have more than one utterance recorded
Extract the MFCC's for all the recordings and store them - these are effectively your "model data" - training samples.

Classification consists of:

Extract the matrix of MFCC's for the test sample
Calculate the DTW of the test sample with each and every MFCC matrix that was extracted and the training stage.
Store the final distance for each comparison. Now you should have a vector of total DTW distances (of a length equal to the number of the training examples):

test vs. class1_sample1: 10.2
test vs. class1_sample2: 20.3
test vs. class1_sample3: 15.1
test vs. class2_sample1: 3.1
test vs. class2_sample2: 2.5
...

For each of the classes (keywords), sort the distances (from small to large)
Pick the k top distances for each class and average them. These are the per-class memberships.
Predict the keyword by choosing the class with smallest distance.


Can the intelligence message of an FM signal be cross-modulated to another center frequency due to non-linearities? Can the intelligence message of an FM signal be cross-modulated to a  a different frequency by non-linearities, but with the FM modulating message remaining (relatively) intact on the new frequency?
This isn't really trying to solve a particular problem, but I'm trying to understand if this is possible, and if the math checks out.  In FM two-way radio, sometimes users experience "bleed over" from another channel.  I've heard intermodulation or cross-modulation being cited as the cause.  Reading the literature however I can't find this phenomenon discussed.  I'm Using Leon Couch, Digital and Analog Communications Systems by the way.
Studying non-linearities; intermodulation and cross modulation, I see that the 3rd order output is:
$$V_{out} = K_{0} + K_{1}V_{in} + K_{2}V_{in}^2 + K_{3}V_{in}^3$$
If you put signals into the system where
$$V_{in} = A_{1}\sin(\omega_1t) + A_{2}\sin(\omega_2t)$$
Then part of the third order outputs are:
\begin{align}3K_{3}A_{1}^2A_{2}\sin(\omega_1t)^2\sin(\omega_2t)  =   \frac{3}{2}&K_{3}A_{1}^2A_{2}\{[\sin(\omega_2t) \\&- \frac{1}{2}\big[\sin(2\omega_1 + \omega_2)t + \sin(2\omega_1 - \omega_2)t\big]\}\end{align}
and another
Generally it's explained that term $\sin(2\omega_1 - \omega_2)t$ accounts for intermodulation distortion that is difficult to filter because it is close to the desired signal.  Also they show that if there is any amplitude variation on signal 1, that it can amplitude modulate signal 2, because
$$\frac{3}{2}K_{3}A_{1}^2A_{2}\sin(\omega_2t)$$
But I never see a discussion of an FM signal being cross-modulated onto another FM signal.  It seems to me that if signal 2 is an FM signal as:   $\sin\big((\omega_2 \pm \Delta\omega) t\big) $. then the cross modulating term:
\begin{align}
\sin(2 \omega_1t - \omega_2t)  &=  \sin\big(2\omega_1t - (\omega_2 \pm \Delta\omega)t\big)\\
&= \sin\bigg(\big(2\omega_1 - (\omega_2 \pm \Delta \omega)\big)t\bigg)\\
& = \sin\bigg(\big(2\omega_1 - \omega_2 \mp \Delta \omega\big)t\bigg)\\
& = \sin\big((2\omega_1 - \omega_2)t \mp \Delta \omega t\big)
\end{align}
Which looks to me like the modulating signal would be impressed on the intermodulation product $\sin((2\omega_1 - \omega_2)t$
It seems to me about the same as super-heterodyning the message signal, which would be the $ \Delta\omega t$
Of course my reasoning might be totally incorrect.
Like I said, I've never been able to find a mathematical or physics treatment of the phenomenon, but I've heard people argue that it happens, and others argue that in the FM domain cross modulation will only cause noise, not the imparting of an intelligent signal on another signal like in AM.

Yes and no,   it can't be intermodded onto another SIGNAL, but it can be intermodded onto another CHANNEL.  For example, if you are monitoring 144 MHz, an FM signal on 145 MHz and on 146 MHz CAN, due to non-linearity, end up being heard on 144 MHz.  But it can't be modulated onto an existing signal on 144 MHz.
Mark


How to calculate EVM in %age of an Equalized Constellation in 16QAM? I have an equalized constellation for 16 QAM. The constellation is equalized by LMS algorithm.

I want to calculate the EVM for the equalized constellation. How can I calculate this so that the EVM lies within the threshold limit?

The Error Vector is the Euclidean distance from the actual sample at the optimum timing location in each symbol to the actual symbol location in a reference constellation (as the distance to closest decision boundary, just prior to decision). The measurement metric Error Vector Magnitude (EVM) is computed as an rms quantity over multiple error vectors, where the waveform and the decision boundaries are scaled to either the rms amplitude or a maximum amplitude in the constellation, depending on which standard is used. We can therefore use one sample for every symbol to compute a sufficient number of errors to derive a statistical quantity.
Below shows a single Error Vector for 16 QAM, and the error vector magnitude would be the root-mean-square of all error vectors as given by the formula above (once the constellation and the received signal are properly normalized as described).  It also must be mentioned that the normalization by the rms of the waveform assumes the symbols in the waveform are equally distributed, but if not normalization can occur through known training sequences. The idea is that the error vector magnitude should be zero if all the symbols land on the exact locations in the constellation.

To compute EVM normalize the constellation (by either the rms, the peak level of an outermost point, or the rms of the magnitudes depending on which definition is used). The waveform samples just prior to decision (one sample per symbol) are similarly normalized which would then minimize the computed EVM (this is consistent with what the receiver would do anyway just prior to decision to minimize error rate, so it is very easy to compute).
Once properly scaled the EVM computation in its general form is just the formula for the standard deviation of samples from a zero mean random process and given as a percentage quantity; so the average error magnitude given as a percentage of a normalized constellation:
$$EVM = \frac{\sqrt{\frac{\sum_n|t_n-\hat t_n|^2}{N}}}{S} \times 100\%$$
(Where $S$ is given by whatever normalization is used for the constellation)
Notice at Wikipedia en.wikipedia.org/wiki/Error_vector_magnitude it is defined as "The ideal signal amplitude reference can either be the maximum ideal signal amplitude of the constellation, or it can be the root mean square (RMS) average amplitude of all possible ideal signal amplitude values in the constellation. ". In this paper they define the normalization as the rms of the constellation: eprints.soton.ac.uk/263112/1/paper_101.pdf , while here the outermost point is used: ieee802.org/16/tg1/contrib/802161c-01_28.pdf .
Point being, if an EVM computation result is being provided or used, it is important to provide the definition and equation that is being used for that particular computation.


Overshoot when using Bessel filter Here is simple code to illustrate overshoot using bessel filter:

sig = (np.random.rand(300)-0.5)*2
b, a = signal.bessel(5, np.asarray([10e3, 400e3])/(0.5*1e6), btype='bandpass')
sig_f = signal.filtfilt(b, a, sig)
print(np.amax(sig_f))
plt.plot(sig, label='input')
plt.plot(sig_f, label='filtered')
plt.legend()

I need a zero phase bandpass filter. I want to reduce overshoot as much as possible, that's why I'm using Bessel filter. But I still have some overshoot in some cases.

In this example, I would like to understand which parameters of my input signal/filter are closely related to overshoot. In my case, is it related to sharp transition in my signal? How can I prevent overshoot?

If it were a low-pass filter, I would say you need to turn the resonance (or $Q$) down enough that both poles are real.

What's the order of your bessel filter?

But a bandpass filter is necessarily resonant.  i might think that a BPF must have poles that are complex.  (now that i think about it, a BPF does not necessarily have poles that are complex conjugate.)

if the BPF has gain exceeding 1 (or 0 dB) for its resonant frequency, then then it should be expected that the output will be larger than the input at that resonant frequency.  but if the gain is less than 1 for all frequencies, i dunno why it would overshoot at that frequency.

but there are zeros for BPF also.  one at DC and another at Nyquist.  so maybe you have to turn down the $Q$ to $\frac12$ so that no pole is "ringing".


On the simplification using trigonometric functions Assume I have a matrix $D$ whose its entries are as below :

Where $A$ and $B$ can be written using using the trigonometric functions for (1) as:

My question, Is it possible to simplify (1)
 more? how can we do that ?

If you use the Euler's formula, you can simplify like this:
$$
[d]_{k,n} = \frac{\sqrt{2}}{N}\left( \cos{\left[ \frac{(k-1)(2n-1)\pi}{2N} \right]} e^{j\frac{2 \pi nk}{N}} \right)
$$

I think we can't simplify more.

PS: If you use your expressions of $A$ and $B$ and again the Euler's formula, you will get the same result.


Neural network with much less positive examples I am developing a medical algorithm working on large amounts of data. However, as in many medical scenarios I have a lot of negative example and only a few positive ones. So any neural network I train is clearly biased towards negative classification. In fact it misclassifies most positive examples and although overall accuracy is OK, false negative rate is sky high.

I wanted to know what is the best way to make the network achieve lower false negative rate, I thought about duplicating the positive examples several times or maybe reducing the negative examples randomly,
What do you think?

There are a couple of things you can do to treat the class imbalance.

A way is to use an SVM-loss (Hinge-loss) function in the network, where the margin around the decision boundary is maximized. The good part of this: Because this function considers the support vectors only, the data points which do not lie close to the margin will not be taken into account. Of course, with more data, you almost guarantee to change the decision boundary and therefore the support vectors. However, especially if your data is spread around the boundary (if you can sample some positive-like true negatives), then you should get close to the correct minimum.
This is an easy way to get around without much effort.
You could augment (not just duplicate) the positives more than you augment the negatives, so that the data becomes balanced.
The weighting of the samples is already discussed, and I think it is valid.


MATLAB gradient derivative troubleshooting I have an array $A$ having the next 145 values

I would like to calculate the $\frac{dA}{dX}$, having a 1D grid, $x$: 1:286:41468

I use the function gradient:

DA_DX = gradient(A, 41468/145)

I am trying to understand why the $\frac{dA}{dX}$ output graph has peaks which cause problems to my code. Please, do you see something wrong in my steps? Is any other way to calculate the $\frac{dA}{dX}$ without gradient function?

You can calculate the derivative using the diff function and then dividing by the time interval between the "samples" in your vector like this:

dA = diff(A) / 286

The diff() function calculates the difference between an element and the previous one in the vector. With this you'll be calculating the derivative for each element in the vector as the increment in A divided by the increment in x. Check out the help for the diff() for more info.


What Does 'Zero Order Hold' and 'First Order Hold' Mean? While studying the Image Magnification in spatial domain, I have come across this definition of Image Magnification by Replication:

  Replication is a zero order hold where each pixel along a scan line is repeated once and then each scan line is repeated.

And the definition for Image Magnification by Interpolation is given as:

  Linear interpolation is a first order hold where a straight line is first fitted in between pixels along a row. Then the pixels along each column are interpolated along a straight line.

Can you please explain the terms 'zero order hold' and 'first order hold' in layman's terms, so that I could understand these two definitions? Thanks.

We need to assume the reader knows some basic stuff to answer that.
Let's give it a try.

Lets understand the sentence - Zero / First Order Hold.
We have the Zero / First Order and the Hold.

Zero / First Order hold means the order of the Taylor Series of the function we use to interpolate. In other words, the degree of the Polynomial we can write the function with.
Zero Order means the function is constant, we interpolate the same value in the missing parts.
First order means we can use linear function to interpolate (Line with a slope).

Hold means we hold the parameters to be the same until the next sample.


What's the difference between "Magnitude Spectrum" and "FFT"? I am particularly talking about the matplotlib's axes.magnitude_spectrum() function that plots the magnitude spectrum. Is this a representation of the discrete Fourier transform?

FFT (fast Fourier transform) is an algorithm that implements the discrete Fourier transform, in a fast manner. Both computes a "full spectrum" from a signal, which is generally complex. The amplitude spectrum is, up to minor corrections (scaling of components), the magnitude of the "full spectrum".


Silence removal from the audio Suppose i've a mp3 file.

I performed FFT of that mp3 file. Then is it good to state the following:
1. Silence occurs by that particular DFT point whose magnitude is lesser
than some threshold. 2. Inverse FFT gives me the original sound file.

1) No, if you have a particularly low level in one of the FFT points it will indicate that you have low energy in that particular frequency, but you can have several other frequencies at the same time.

2) In theory yes, but in practice the numerical precision will make the resulting signal slightly different than the original one.

Note: before applying the FFT you'll need to decode the mp3 file in order to get the audio data.