Compress a signal by storing signal diff instead of actual samples - is there such a thing? I am working with EMG signals sampled at 2kHz and 16 bits, and noticed that they "look smooth", that is, the signals are differentiable, and if I apply a "diff" function (numpy.diff in my case) the magnitude of the values is considerably lower than the actual samples.

So I am considering to do something like:

Split the signal into chunks of a given size;
Foreach chunk, using variable length quantity (or similar), create a byte list and:

For the first sample of the chunk, add its absolute value;
For the remaining samples of the chunk, add their difference, relative to the previous value;

This way, the smoother the signal, and the closer it is to the baseline, the more I expect to decrease the byte-size of each chunk, by decreasing the individual byte-size of a large part of the samples.

Although I suspect this would improve things for me, I also suspect that this is nothing new, and perhaps it has a proper name, and even more elegant/efficient ways to implement it.

So the question is: what is the name of this compression technique, and what are its alternatives and/or variants?

You can also think of delta encoding as linear predictive coding (LPC) where only the prediction residual ($x[n]-\hat{x}[n]$ in @robertbristow-johnson's notation) is stored and the predictor of the current sample is the previous sample. This is a fixed linear predictor (not with arbitrary coefficients optimized to data) that can exactly predict constant signals. Run the same linear predictive coding again on the residuals, and you have exactly predicted linear signals. Next round, quadratic signals. Or run a higher-order fixed predictor once to do the same.

Such fixed predictors are listed in Tony Robinson's SHORTEN technical report, yours in Eq. 4, and are also included in the FLAC lossless audio codec although not often used. Calculating the best prediction coefficients for each data block and storing them in a header of the compressed block results in better compression than the use of fixed predictors.

For $m$-bit input the residual is an $m+1$ -bit number, because it is the difference of an $m$-bit input and an $m$-bit prediction. However, removing the most significant bit (MSB) of the residual has no consequence in $m$-bit modular arithmetic, so the residuals can be stored as $m$-bit numbers.

The linear predictor is supposed to do the whitening, making the residuals independent. In lossless compression, what is left to do is to entropy code the residuals, instead of using run-length or other symbol-based encoding that doesn't work so well on noisy signals. Typically, entropy coding is done by a prefix code (also known as prefix-free code) that assigns longer code words to large residuals, approximately minimizing the mean encoding length for an assumed distribution of the residual values. A Rice code (also known as Golomb–Rice code or GR code) variant compatible with signed numbers can be used, as is done in FLAC (Table 1), or signed exp-Golomb code as is done in the h.264 video compression standard. Rice code has a distribution parameter that needs to be optimized for the data block and saved in the block header.

Table 1. Binary codewords of 4-bit signed integers encoded in Rice code with different Rice code parameter $p$ values, using FLAC__bitwriter_write_rice_signed (source code). This variant of Rice code is a bit wasteful in the sense that not all binary strings are recognized as a codeword.
$\begin{array}{rl}
\begin{array}{r}\\-8\\-7\\-6\\-5\\-4\\-3\\-2\\-1\\0\\1\\2\\3\\4\\5\\6\\7\end{array}&\begin{array}{lllll}
p=0&p=1&p=2&p=3\\
000000000000001&000000010&000110&01110\\
0000000000001&00000010&000100&01100\\
00000000001&0000010&00110&01010\\
000000001&000010&00100&01000\\
0000001&00010&0110&1110\\
00001&0010&0100&1100\\
001&010&110&1010\\
1&10&100&1000\\
01&11&101&1001\\
0001&011&111&1011\\
000001&0011&0101&1101\\
00000001&00011&0111&1111\\
0000000001&000011&00101&01001\\
000000000001&0000011&00111&01011\\
00000000000001&00000011&000101&01101\\
0000000000000001&000000011&000111&01111\end{array}\end{array}$

As a further enhancement, encoding not just one but multiple residuals into a single codeword can more accurately accommodate the true distribution of residuals and may give a better compression ratio, see asymmetric numeral systems.


At what time stretching ratio would an audiophile notice a wrong pitch I have some low quality audio from my native tongue on a VHS. I am using this audio to dub HQ released of the same program (from another language source). The video frames/second ratio is 25:23.98 . I do not want to tamper with the fps ratio, since I am still using the other audio tracks as alternative languages in my video.

I can't notice a pitch difference from this conversion, I can only hear the pitch-changing VHS artefacts.

Is a 25/23.98 ratio enough to need pitch correction? What is the rule of thumb? Would the VHS artefacts hide the pitch difference (too much pitch-noise to actually matter) or would the average pitch difference drive an audiophile up the wall?

Audiophiles don't hear anything unless they are told what to expect.

Joke aside, this is equivalent to transposition by a bit less than a semitone. The effect won't be noticed on voices - this is way too low for any kind "chipmunk" effect to be observed. As for music, some people very familiar with the original music might detect the change of tonality and/or tempo. But such stretching ratios are not unusual at all - DJs doing beat-matching, TV channels speeding up a sitcom to squeeze in more ads, and many, many movie to video conversions (telecine) have routinely used it.


What is an example of using aliasing to your advantage when recovering an input signal? Suppose you have an arbitrary analog input signal $x_a(t)$ guaranteed to have frequencies within a bandwidth $[f_1,f_2]$ Hz.

Suppose your sampling frequency $F_s$Hz, and sample $x_a(t)$ to produce $x(n)=x_a(\frac{n}{F_s})$.

Then by the sampling theorem, you can successfully recover (or reconstruct) the signal so long as $f_2 < \frac{F_s}{2}$. Otherwise you will experience aliasing and cannot recover the original signal $x_a(t)$.

It appears (from this article) that there are "tricks" that can be employed that take advantage of aliasing scenario and can still recover the signal $x_a(t)$ completely.

The specific part of the article of concern is pasted below:

  Using Nyquist aliasing as benefit

  The trick is to use the aliasing (or frequency folding) to your
  advantage. By undersampling the data converter, higher-frequency
  content will be aliased into all of the lower Nyquist zones (see
  Figure 2). You will need to make absolutely sure that nothing ends up
  in the lower bands – any noise or frequency components in the lower
  zones will also be aliased into the first Nyquist. The good news is
  that the data rate from the data converter is only a fraction of the
  required RF input sample rate if this were a first Nyquist system.
  Under sampling greatly reduces the data rate of the samples supplied
  to the digital signal processor (DSP) or FPGA.

Could somebody explain this to me? Why is the Texas Instruments link not violating Nyquist's Theorem?

you can successfully recover (or reconstruct) the signal so long as $f_2<F_s/2$

Nope. The sampling theorem says that you need 2 samples per Hz of bandwidth, so in this case you'd need

$$Fs > 2*(f_2-f_1)$$

For more info on how this works google "bandpass sampling" or just read https://en.wikipedia.org/wiki/Undersampling.

The basic idea is that sampling creates a periodic repetition of the original spectrum with a period of $F_s$. As long as none of the repeated spectra falls on the same frequencies as the original spectrum, you don't get aliasing.


how to generate a smooth sine wave with time varying frequency in matlab I'm currently designing a knock sensor for my FSAE team's car and I'd like to be able to test various software strategies with a fairly realistic signal (at least more realistic than a couple sine waves added together with some white noise).

From some reading, the frequencies I'll be getting will all be related to a fundamental frequency directly related to engine speed. As engine speed varies continuously, I'd like to be able to take a vector of times and transform it into a vector of frequencies that smoothly and randomly transition within the range 20-170 Hz, ideally with a uniform distribution. I can then create a couple sine waves based on these frequencies and a couple multiples, add them, add some white noise, add a knock at random points, then feed that into my knock detection code.

It's generating a smooth, uniformly distributed frequency vector that is my current issue, but I'm also open to other solutions to the problem.

The wording of your question (your terminology) makes it difficult to understand exactly what you desire. Is your question, "Using Matlab, how do I generate a time-domain sequence representing a sinusoidal wave whose instantaneous frequency increases as a linear function of time so that its spectrum contains a broad band of spectral energy?" If so, then try the following example Matlab code:

Fs = 1000;       % Signal sample rate in Hz

Dur = 1;         % Signal's 1-second time duration

Time = 0:1/Fs:Dur;

Fstart = 100;    % Start freq is 100 Hz

Fstop = 400;     % Final freq is 400 Hz

Sig = chirp(Time,Fstart,Dur,Fstop);   % Generate a "chirp" sequence

rpatel, generate the above 'Sig' signal, compute its DFT samples, and plot the DFT's spectral magnitude values. Maybe that's what you want.


Speech Overlap Detection I’m trying to figure out how to go about differentiating between a single-voice speech and a multiple-voice conversation in which voices overlap. In other words – the moment of overlap is what I’m interested in. I want to be able to listen to a vocal conversation between two people, and detect when voices collide/overlap.

I guess that spectrum-wise there’s a major difference between a speech-based conversation with no overlaps compared with a conversation with speech overlaps, so perhaps one direction to look at is detecting sudden spectral changes?

Most human vowels produces a series of harmonic spectral peaks that resolve to a single F0 pitch (plus some broad spectrum noise).  If you can't fit a single F0 pitch to the series of the spectral peaks, then you might try testing the hypothesis that the spectrum resolves to 2 or more f0 pitches, which would imply more than one speaker talking with an overlap of vowel sounds in their words.

Detecting an overlap between a vowel and someone else's consonant, without lots of false positives, seems like it might be far more difficult.


Find Magnitude Responce at $\omega=\dfrac{\pi}{2M}$ where $M$ is the number of channel in QMF I have frequency response of low pass filter

    M=8  %%%channel
    x=8*512
    d1=fir1(N-1,wc/pi,w);
    % L=length(d1);
    % for k=1:L
    %     b2(k)=((-1)^k)*d1(k);
    % end
    %%%%%%Frequency response
    h1=fft(d1,x);
    a1=max(abs(h1));
    m1=20*log10(abs(h1/ a1));
    c=0:x-1
    plot(f1/x,m1)

I have to find magnitude responce at $\omega=\dfrac{\pi}{2M}$ i.e., $h_1$ at $\omega=\dfrac{\pi}{2M}$

Your code has too many unassigned variables for us to execute it. Your FFT freq axis, in radians/sample, goes from 0 –to- 2pi. The equivalent mathematical FFT integer freq axis index goes from 0 –to- 4095. You want to know h1 at w = pi/(2M) radians/sample.  Let's use variable k to represent the mathematical FFT integer index corresponding to your w.

So we can set two ratios equal to each other: k/4096 = w/2pi.  Next we write k/4096 = (pi/(2M))/2pi. This gives us k = 4096/(4M) = 4096/32 = 128 as the math index corresponding to your w. Because of Matlab's unpleasant indexing method the magnitude response you want to examine is your h1(k+1) = h1(129) FFT sample.  [By the way, as long as your M is an integer power of two, k will be an integer. If M is not an integer power of two, you'll have to use either a 'floor()' or 'ceil()' command to compute the mathematical FFT integer index k.]


Does a PID need to be "online" to be tested? Can the output signal be viewed without it being used? I am taking my first shot at creating and tuning a PID for my software application. The PID I am creating is used to control the temperature of a "thing". The thing is heated by focusing a laser on it and changing the laser's output power in the range of ~0.1W to 1.0W. The temperature is then calculated and is in the range of ~100C to 1000C.

An example sentence/problem to refer to. I desire the temperature of my thing to be 600C, to start I set the laser's power to 0.5W, then measure a temperature of 500C, the PID is iterated and calculates a output signal of 0.6W (Assuming the PID parameters are tuned correctly).

// Is this the correct terminology?
???            = Laser power of 0.5W
Setpoint       = Desired temperature of 600C
Measured value = Measured temperature of 500C
Output         = Output signal of 0.6W

To test my PID implementation and parameters I have been allowing the PID to iterate and calculate an ouput signal, but I have NOT been allowing this signal to command the system ( In the example I have been calculating the output signal of 0.6W but not setting the laser power to 0.6W; it stays at a constant 0.5W ). This has somewhat worked to hone in on appropriate PID parameters but I think I need to put the PID "online".

Obviously my PID parameters, and maybe implementation, are not going to be correct until they are refined and tuned so I do not want to allow it to control the system for various types of concerns and safety issues it could cause. Is this ok to do? Does it make sense?

I think this is a flawed way to see if the output signal is correct because the Error and Integral from past iterations are going to reflect past commanded signals. The error/integral could continue to grow even though they shouldn't be?

EDIT: Rephrased. Will a PID always produce a sane output whether it is online or not?

The PID can be tested in an open-loop setup, however any integral gain <> 0 will cause the PID output to run away, so you are likely to hit clamping (like control voltage maxima) with your test rig.

What really matters and needs to be tuned is the closed-loop behaviour. This includes your plants behaviour. Stability is prime goal, maximum overshoot and settle time (to set-point change or disturbance) can be tuned.

If you cannot tune in the live system, modeling and simulation is the way to go (look for Simulink Control Design). You would typically start with an empirical plant model, tune your PID for it, then run the actual device with this PID but a set-point well below Tmax. Instrumentation will give you data to improve your plant model. Repeat this development cycle a few times until you are confident to increase to the actual set-point.


What is the absolute phase and the relative phase, of a signal? I need to know what is the absolute phase and the relative phase, of a signal?
and why these phases are important in the signal processing?

"Absolute phase" is really relative phase to a point in time at which the phase measurement is referenced (Noon GMT on Tuesday, the first sample in a time-domain window of sample data for instance, or maybe the middle sample in a vector).
Relative phase can also be measured with respect to another waveform of the exact same constant frequency, as the difference in phase measurement between the two waveforms at any one point will be the same at all other points, as long as the two frequencies stay identical.

Phase is important because it is one of the 3 parameters needed in an equation that completely describes a sinusoid.  Without phase information, the waveform is ambiguous.

$$y(t) = a\cdot \cos\left(b\cdot t + \phi_{t_0}\right), \quad \text{ where }\phi_{t_0} \text{is the phase at time } t = t_0.$$


Understanding the mathematical proof for the alias frequencies in a sampled sine wave I'm struggling to get my head round the mathematical proof for the alias frequencies in a sampled sine wave.

I understand that sampling a sine wave of frequency $f_0$ every $t_s$ seconds gives you:

$$x[n]=\sin(2\pi f_0nt_s)$$

I also understand that, because the sine wave is periodic every $2\pi$, you can add any multiple of $2\pi$ to the angle and get the same values for the sine, i.e.,

$$\sin(2\pi f_0nt_s)=\sin(2\pi f_0nt_s+2\pi m) \quad\text{(where $m$ is any integer).}$$

The proof I'm looking at then factors out $2\pi$ and $nt_s$ to get:

$$\sin\left(2\pi(f_0+\frac{m}{nt_s})nt_s\right)$$

...but then it says to let $m$ be an integer multiple of $n$ so we can replace the $\frac{m}{n}$ ratio with an integer $k$.

I don't understand how $m$ can go from being "any integer" to "an integer multiple of $n$". If $m$ is any integer and $n$ is an integer then how can the ratio between them be an integer?

I know I'm missing something obvious here and I'm searching for that light-bulb moment but it's not happening. Because this is so fundamental to DSP I don't just want to accept the formula and move on without thoroughly understanding it.

The reason is that if it is true for any $m$, it is also true for $m=kn$.

I will sketch the proof in another way.

Call $f_s = 1/t_s$ sampling frequency where $t_s$ is sampling period, the two signal $x(t) = \sin(2\pi f_0 t)$ and $x_k(t) = \sin(2\pi (f_0 + k f_s) t)$ have the same values at sampling instants (aliasing), i.e. $x[n] = x_k[n]$.

Indeed,

\begin{align}
x[n] &= \sin (2 \pi f_0 n t_s) \\
x_k[n] &= \sin \left(2 \pi (f_0 + k f_s) n t_s\right) \\
&= \sin (2 \pi f_0 n t_s + 2\pi k f_s n t_s) \\
&= \sin (2 \pi f_0 n t_s + 2\pi k n) = \sin (2 \pi f_0 n t_s)\\
&=x[n]
\end{align}


Filter Design - Step Response I'm  looking for a tool\method which allows me to design a filter without involving the frequency response at all.
Specifically, I'm trying to design a simple FIR smoothing filter for an emmbedded system which reads different analog sensors. My sampling rate is 1kHz and I'd like to generate an FIR filter with 0-overshoot and a smoothing effect (e.g. in attached pic).

The different tools I found (including MATLAB) all involved a design around a LP filter, which means I have to analyze the frequency response of my system.

Any relevant tools\methods out there?

For linear filters, binomial filters, which can be considered as FIR approximation to Gaussians, are very simple. Their coefficients are given by:
$$ h_k=\frac{ \binom{n}{k}}{2^n} \,.$$
For instance, you can get $H_3 = [ \frac{1}{4},\frac{1}{2},\frac{1}{4}]$, or $H_5 = [ \frac{1}{16},\frac{1}{4},\frac{3}{8},\frac{1}{4},\frac{1}{16}]$. Since their coefficients are positive, they do not overshoot.

Here are graphs for $n = 5,17,25,31$:

Savitzky-Golay filters, performing polynomial interpolation, are interesting too, but tend to   overshoot with higher degrees.

Do not fool yourself though: as long as the filters are linear, there is a notion of frequency very close behind.


modulate/demodulate audio file as FSK signal mixed with music signal What I want to do:
I want to encode an audio message inside a regular music file; this shall be decoded by an audio analyzer running on another machine (by listening to the audio playback from the first machine).

What I have tried so far:
Encode:

Change the pitch of secret audio message to D#/EB-10 (which is 19.5
KHz according to this site:
http://www.flutopedia.com/pitch_to_frequency.htm). This should
convert the secret audio message to something that falls outside of
the human audible range.
Merge this secret audio message with music file.

Decode:

Apply high pass filter at 19 KHz.
Change pitch back to the original range (e.g. B-6)

I'm using Audacity to do the above steps.

The problems that I'm facing:
1. The "secret" message is still audible in the encoded output.
2. Decoded output is significantly different than the secret message.

I'm new to the audio dsp domain, so I'd appreciate a review of my approach and suggestions to achieve the goal that I described above.

Edit:
Below is the spectrum plot for "secret" audio message generated using Audacity.

Edit2:
Below image shows linear spectral chart for "secret" signal mixed with generated tone of 19kHz.

Edit3:
The below spectral graph is for the carrier "music" file.

So, your unmixed secret signal has energy content in a frequency range of ca -6 kHz to +6 kHz.

As far as I understand, you mix that with a real cosine of ca 19 kHz, so that it now has signal content in -25 kHz to -13 kHz, and +13 kHz to +25 kHz. All things are still nice and symmetric – so if you mix it back down to 0 Hz, you'd still have the original spectrum.

Only that your 48 kHz sampling rate doesn't even allow your to unambigously represent a signal above 24 kHz – it instantly aliases into the lower end, and vice versa, so you have distortion in your signal of interest (still shouldn't be hearable).

But then you high-pass filter at 500 kHz above the frequency you shifted to – but that cuts off 5.5 kHz of your signal! So, afterwards, your signal is no longer symmetrical in spectrum, and mixing it back down will produce non-constructive aliases.

So, to understand what's happening here, plot your secret signal's spectrum after being shifted up. Don't use the logarithmic frequency axis, but simply linear frequency.


Finding values of DTFT without explicitly computing My attempt :

a) Summation of all values?

b)c)d) Failed

e) Parserval's theorem

You just need to know the formulas for the DTFT and its inverse:

$$X(e^{j\omega})=\sum_{n=-\infty}^{\infty}x[n]e^{-jn\omega}\tag{1}$$

and

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

From $(1)$ you see that your answer for $(a)$ is correct. Also the answer for $(d)$ follows immediately, if you realize that $e^{-jn\pi}=(-1)^n$.

For $(b)$ it's important to see that the coefficients are asymmetrical. Looking at $(1)$ you see that you can pair the terms for indices $\pm 1$, $\pm 2$, etc., which allows you to rewrite the sum $(1)$ as a sum of weighted sines time $j$ (imaginary unit), because $e^{jx}-e^{-jx}=2j\sin(x)$. So you'll have a purely imaginary expression. Computing the phase should then be easy.

$(c)$ is easily solved using formula $(2)$ (which value of $n$ do you need to plug into $(2)$?). And finally, for $(e)$ you're right that you need Parseval's theorem.


What is the algorithm to generate sine waves of arbitrary frequency in the STFT domain? I'd like to write a DSP algorithm to do additive synthesis using arbitrary sine waves with inverse rectangular FFTs.
This requires two things:

The ability to generate phase/amplitude lists that cause the IFFT to create sine waves that are not integer multiples of 1/(IFFT length)

Example:  How do I generate 1000 Hz with 16 sample IFFTs at 44100 Hz?

The ability combine two frames pre-IFFT so they are added post-fft.

What are the necessary steps here?

The FFT of a non-integer multiple frequency is very complex, you don't want to go there!
There are various oscillator schemes that are a few MAC instructions per sample such as https://vicanek.de/articles/QuadOsc.pdf
Checkout a review here: https://www.njohnson.co.uk/pdf/drdes/Chap7.pdf


Causal Noise Free 1st Order Derivative in Discrete Domain I need to have causal noise free first order derivation (Derivative).
Now I am using a simple finite differences formula:

$$ \frac{ x \left( n \right) - x \left( n - 1 \right) }{ {T}_{s} } $$

The filter above, being High Pass, is super sensitive to noise.
Do we have a noise free or noise reduced first order differentiation which does not need more than three point for calculation?

Well, Since basically the Derivative Operation is Linear Filter you can chose your own optimal trade off between Noise Sensitivity and Bandwidth.

If you look at Finite Differences Coefficients page at Wikipedia you can see you can chose higher "Accuracy" filters for 1st Derivative.
You can also chose Forward / Central / Backward method.


Should i high pass audio for real time frequency detection? Haven't been able to find a clear answer on this. I'm interested in 75-500Hz bandwidth. It's for a tuner(android app) so i want a pretty quick response time.

Is it worth the computational time to high pass the signal?

Yes, it is. If you pickup signals with a microphone you will often see enormous amounts of low frequency signals (air condition, HVAC, traffic, wind, etc.). Even in quiet rooms there tends to be a lot of very low frequency noise.

This noise can dominate the time domain waveform and reducing it can often make your pitch detector more reliable and/or allows for a cheaper or more efficient algorithm.


Camera calibration using multiple images I am trying to get the hang of camera calibration. For this purpose, I have mainly watched this lecture on youtube: Lecture 12: Camera Model by UCF CRCV. I believe that I understood the main ideas:

We basically deal with the question to relate 3D world points to 2d image points. In order to do this

We take the 3D point W in homogenous coordinates (in a defined world coordinate system)
Apply transformations (rotation, translation) to bring the camera and world coordinate system in alignment
Use perspective projection (from the pinhole camera model) to calculate the image point I.

If we do not have the perspective projection (that is the intrinsic parameters) and we do not have the camera position and orientation in the world (extrinsic parameters), we can calculate them (=calibration) by using an object with easily findable points in the world, identifying those points in the taken image and building a system of equations which we can then solve using least squares.

I now noticed, that many calibration tools actually require multiple images to calibrate reliably. Concerning this, I have a couple of questions:

How do we relate the multiple images we might have taken from the same object that we have moved? I have seen that the world origin is often placed in a chessboards corner, but this origin will change if you move the chessboard. Don't the origins have to match, such that we can calculate the extrinsic parameters over multiple images?
Why exactly is it better to use multiple images? Couldn't you just use a big chessboard that gives you a sufficient number of points to solve the system of equations?

You are correct: to calibrate a camera you need a correspondence between 3D world points and 2D image points. The problem is that the 3D points cannot be co-planar, so people were building 3D calibration rigs, e. g. a box made of checkerboards. One image of a rig like that would be enough to calibrate, but those rigs are hard to build, because you have to get the planes to be at exactly right angles to each other.

Then Zhengyou Zhang came up with the algorithm to calibrate from multiple images of a planar pattern, like a checkeboard. You still need non-coplanar points, but you get them from multiple images of a planar pattern in different 3D orientations, rather than a single image of a 3D rig.

Noise and errors are another reason why you need multiple images. You cannot detect the checkerboard corners with infinite precision. So you need a lot of points, to estimate the camera parameters more reliably.


Removing empty/low amplitude part from signal I have a signal of length 'N' which looks like this

[.......p-set...] [SINE-WAVE] [....q-set...]

amplitude of points in p-set << amplitude of points in sine-wave

Sometimes there are a few sparse points in p-set whose amplitude is comparable to sine-wave points. However, number of such points is much less than the length of the sine wave.
I want to find the beginning point of the sine-wave.

Clearly I cannot use thersholding as those sparse points will be counted as the beginning of the sine-wave which is wrong.

Is there any effective solution or existing functions? Sliding window may be? Please explain with some code (MATLAB). Due to firewall restrictions I cannot upload image. Sorry.

If you know the frequency of the sine signal  already, you can correlate your [.......p-set...] [SINE-WAVE] [....q-set...] signal with the sine signal. The correlation will give you a peak where sine wave exactly matches with the signal. Since you know the sine signal length you can figure out the beginning of the sine signal in your original signal. I hope this helps you a bit.


Shifting Sequences Using MATLAB I'm reading the book
Vinay Ingle, John Proakis - Digital Signal Processing Using MATLAB® (3rd Edition) and I came across a function to shift sequences. Mathematically,

$$
y(n) = \{x(n-k) \}
$$

if $m = n-k$, then
$$
y(m+k) = \{ x(m) \}
$$

function [y,n] = sigshift(x,m,k)
n = m + k;
y = x;
end

This function shifts nothing. I've tried to find errata list of the book, but I couldn't find it. Is this function correct?

In order to shift in MATLAB you need to play with the indices of the signal vector.
for example, given vector x - vX of length L, to shift it in time such that the first sample is the fifth you'll do this:

vXShifted = x(5:L);

Pay attention that the signal is shorter by 4 samples.

Now, sometimes, for periodic signals, a cyclic shifting is required which can be done using circshift() function in MATLAB.


How can I generate a pulse train in the frequency domain? I just need to generate the equivalent train in the frequency domain by hand.

I know that the Fourier Transform of a pulse train is a pulse train, with the intervals of the pulses changed by (1/T).

However I'm confused on the implementation when you are also working with sampling rate of the original signal.

Can anyone quick point me in the correct direction.

Thanks!

The train of pulses in the time domain:
$$x(t)=|||(t)=\sum_{i=1}^{\infty} \delta(t-i)$$

Equals to the following train of pulses in the frequency domain:
$$\mathbb{F}x=|||(f)=\sum_{i=1}^{\infty} \delta(f-i)$$

Where $|||$ is the train of unitary pulses, comb standard signal.

Hence:
$$|||(t) \rightarrow |||(f)$$

The train of unitary pulses in time domain separated by $\tau$ is (do you know why we are having $1/\tau$?):
$$1/\tau|||(t/\tau) \rightarrow |||(\tau f)$$

And the train of unitary pulses in frequency domain separated by $1/\tau$ is:
$$|||(t/\tau) \rightarrow \tau|||(\tau f)$$

NOW when having a sampling time of $dt$, each delta pulse in the instant $i$ can be represented in discrete time as:
$$\delta(t-i)=\frac{1}{dt}\delta_{ti}$$
Where $\delta_{ti}=1$ for $t=i$ and 0 elsewhere, as in the discrete kronecker delta sense, approximating each pulse as the smallest possible rectangle (other approximations are possible though...).

In discrete time/frequency, we just have to replace the proper $delta$ implementation, as already defined.

Finally, the $\infty$ in the time sum is bounded to $T$, the total measuring time, and to $F=0.5/dt$, the sampling frequency band (nyquist?). So the period $\tau$ must be divisor of $T$, $T=\phi\tau$, and $1/\tau$ must be a divisor of $F=0.5/dt$, $T\phi=0.5\tau N$  for the numbers in the FFT algorithm to be squared in order to create a perfect periodic sequence in both the time and frequency streams.

HENCE, the standard FFT algorithm, for example with $\tau=6$, $\phi=3$ and $dt=1$ (for simplicity):

tau=6;
phi=3;
T=phi*tau;
x=1/phi*repmat(eye(1,tau),1,phi);
Fx=fft(x)

Implements the non unitary train of pulses in time domain:

x=1/3*[1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0]

with the unitary train of pulses in frequency domain:

Fx=[1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0]


Sine sweep vs impulse response What kind of pros and cons does these method for identifying a system have?
I seem to have a hard time finding literature discussion in this issue.

I am operating on a linear system, or the linear area of the system. So non linearity isn't an issue here.

A good starting point to better help you understand the differences between several available methods for Transfer Function measurements is the Müller and Massarani paper entitled "Transfer-Function Measurements with Sweeps".  Whilst the title eludes to coverage of only the ESS method, it does do a good job at covering other techniques too.

Angelo Farina has also contributed many works regarding Transfer Function measurements for a variety of applications, all of which he makes freely available here:  http://pcfarina.eng.unipr.it/Public/Papers/list_pub.htm

I hope that this assists with your literature quest.  Following reading these papers you should have a better understanding of the different methods, along with their respective pro's and con's, so you can make an educated decision regarding which technique best suits your application.


Root raised cosine parameter I try to create matched filter for receive BPSK signal. Sampling rate of the signal is 228000 Hz, data rate 1187.5 symbols/sec, signal was shaped using root raised cosine filter with rolloff factor = 1.

I cant get what should be sps — Samples per symbol in function rcosdesign . First i thought it's 228000/1187.5 , but tried it and probably I'm wrong...

No, you're right. It wants to know how many samples you have for each symbol in the BPSK stream. This is the ratio of the sample rate to the symbol rate.


Offset QPSK detection in GNU Radio - Sample delay An OQPSK detector is being tested in GNU Radio. The architecture was obtained from Michael Rice's Digital Communications - A discrete-time approach. The flowgraph is shown below.

The modulator architecture is basically QPSK delayed by half the sample rate, as shown below. In the flowgraph, the block oqpskIQMap maps the I/Q symbols from QPSK to OQPSK. It does so by delaying the Q-channel by d_delay, which is equal to the number of samples per symbol divided by two.

oqpskIQMap_impl::oqpskIQMap_impl(int sps)
      : gr::sync_block("oqpskIQMap",
              gr::io_signature::make(1, 1, sizeof(gr_complex)),
              gr::io_signature::make(1, 1, sizeof(gr_complex))),
        d_delay(sps/2)
    {
        //set_history(d_delay);
    }

    /*
     * Our virtual destructor.
     */
    oqpskIQMap_impl::~oqpskIQMap_impl()
    {
    }

    int
    oqpskIQMap_impl::work(int noutput_items,
        gr_vector_const_void_star &input_items,
        gr_vector_void_star &output_items)
    {
      const gr_complex *in = (const gr_complex *) input_items[0];
      gr_complex *out = (gr_complex *) output_items[0];
      for(int i = 0; i < noutput_items; i++)
          out[i] = gr_complex(real(in[i]),imag(in[i - d_delay]));

      // Do <+signal processing+>

      // Tell runtime system how many output items we produced.
      return noutput_items;
    }

The receiver architecture from Rice's book is shown below. In this architecture, the match-filter produces 2 samples per symbol. The oqpskIQDemap block processes the output samples (X(KTs), X(KTs + Ts/2), Y(KTs) and Y(KTs + Ts/2) into a constellation point (X(KTs) ,Y(KTs + Ts/2)), effectively reducing the sample rate to 1 sample per symbol.

oqpskIQDemap_impl::oqpskIQDemap_impl(int delay)
      : gr::sync_decimator("oqpskIQDemap",
              gr::io_signature::make(1, 1, sizeof(gr_complex)),
              gr::io_signature::make(1, 1, sizeof(gr_complex)), 2),
        d_delay(delay)
    {}

    /*
     * Our virtual destructor.
     */
    oqpskIQDemap_impl::~oqpskIQDemap_impl()
    {
    }

    int
    oqpskIQDemap_impl::work(int noutput_items,
        gr_vector_const_void_star &input_items,
        gr_vector_void_star &output_items)
    {
      const gr_complex *in = (const gr_complex *) input_items[0];
      gr_complex *out = (gr_complex *) output_items[0];

      for(int i = 0; i < noutput_items; i+=2)
          out[i] = gr_complex(real(in[i]),imag(in[i + 1]));

      // Tell runtime system how many output items we produced.
      return noutput_items;
    }

The transmitted constellation (top) looks pretty okay. The Rx constellation (bottom), on the other hand, seems to have some points crossing the boundaries unexpectedly, given that the SNR is 25 dBs. I suspect the problem is with the way I introduce the delay in the two IQ map/demap blocks.

Please do let me know what you think.

Regards,

Thanks to Marcus' answer, I was able to resolve the problem. See below for the updated code:

IQ Mapping

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <gnuradio/io_signature.h>
#include "oqpskIQMap_impl.h"

namespace gr {
  namespace oqpsk {

    oqpskIQMap::sptr
    oqpskIQMap::make(int sps)
    {
      return gnuradio::get_initial_sptr
        (new oqpskIQMap_impl(sps));
    }

    /*
     * The private constructor
     */
    oqpskIQMap_impl::oqpskIQMap_impl(int sps)
      : gr::sync_block("oqpskIQMap",
              gr::io_signature::make(1, 1, sizeof(gr_complex)),
              gr::io_signature::make(1, 1, sizeof(gr_complex))),
        d_delay(sps/2)
    {
        //set_history(1);
    }

    /*
     * Our virtual destructor.
     */
    oqpskIQMap_impl::~oqpskIQMap_impl()
    {
    }

    int
    oqpskIQMap_impl::work(int noutput_items,
        gr_vector_const_void_star &input_items,
        gr_vector_void_star &output_items)
    {
      const gr_complex *in = (const gr_complex *) input_items[0];
      gr_complex *out = (gr_complex *) output_items[0];

      for(int i = d_delay; i < noutput_items; i++)
          out[i] = gr_complex(real(in[i]),imag(in[i - d_delay]));

      // Do <+signal processing+>

      // Tell runtime system how many output items we produced.
      return noutput_items - d_delay;
    }

  } /* namespace oqpsk */
} /* namespace gr */

IQ Demapping
Please note that this is a downsampling block and which expects 2 samples per symbol and produces 1 sample per symbol (a single constellation point).
    #ifdef HAVE_CONFIG_H
    #include "config.h"
    #endif

#include <gnuradio/io_signature.h>
#include "oqpskIQDemap_impl.h"

namespace gr {
  namespace oqpsk {

    oqpskIQDemap::sptr
    oqpskIQDemap::make(int delay)
    {
      return gnuradio::get_initial_sptr
        (new oqpskIQDemap_impl(delay));
    }

    /*
     * The private constructor
     */
    oqpskIQDemap_impl::oqpskIQDemap_impl(int delay)
      : gr::sync_decimator("oqpskIQDemap",
              gr::io_signature::make(1, 1, sizeof(gr_complex)),
              gr::io_signature::make(1, 1, sizeof(gr_complex)), 2),
        d_delay(delay)
    {
        set_output_multiple(2);
    }

    /*
     * Our virtual destructor.
     */
    oqpskIQDemap_impl::~oqpskIQDemap_impl()
    {
    }

    int
    oqpskIQDemap_impl::work(int noutput_items,
        gr_vector_const_void_star &input_items,
        gr_vector_void_star &output_items)
    {
      const gr_complex *in = (const gr_complex *) input_items[0];
      gr_complex *out = (gr_complex *) output_items[0];

      //for(int i = 0; i < noutput_items; i+=2)
          //out[i] = gr_complex(real(in[i]),imag(in[i + 1]));
      out[0] = gr_complex(real(in[0]),imag(in[1]));

      // Tell runtime system how many output items we produced.
      return 1;
    }

  } /* namespace oqpsk */
} /* namespace gr */

Regards,


Introduction to Statistical Signal Processing There is a postdoc in my lab whose specialty is "statistical signal processing". He has a PhD in Electrical Engineering and he analyzes the neural data collected.

I am wondering what courses/topics I should start studying to follow in his footsteps. I'm not exactly looking for things like statistics and signal processing, I've had basic classes in both but still find it hard to understand his work.

Sometimes there are courses entitled 'statistical signal processing', that's a good place to start :-)   If your university doesn't have this, try looking for 'detection and estimation', or 'advanced signal processing'.  If you don't have a university handy, you could try http://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-432-stochastic-processes-detection-and-estimation-spring-2004/

Much statistical signal processing is linear, so you should learn as much linear algebra as you can. Stocastic processes is a foundational course.  Control theory shares much with SSP, and would be very useful.

This should be enough for a start :-)


Unit Step Response of an Integrator Desinging a very basic 1st order integrator with the continous transfer function:

$H(s) = -1 / (0.001 s)$
[No zeros, one pole at origin]

I was expecting the unit step response to be positive-slope.
However in my Python code I get a negative slope.

num = (0 , -1)
den = (0.001 , 0)
Hs = sig.lti(num,den)
t, s = sig.step(Hs)
plt.plot(t,s)

Provides:

Your integrator has a negative sign:

$$H(s)=-\frac{1000}{s}\tag{1}$$

That's why its step response also has a negative sign:

$$y(t)=-1000\int_{-\infty}^tu(\tau)d\tau=-1000\cdot t\cdot u(t)\tag{2}$$

where $u(t)$ is the unit step.


Any references for the trade-off between software and mechanics/optics in industrial machine vision systems? I explain my question with a simplified example.

I can design an industrial vision system for the automatic inspection of an item with these main requirements:

the image of a good piece must be a black background and the piece must be grey.
the defect must appear as a white area inside the grey area.

These requirements simplify a lot the software part of the system: in order to classify an item as defective the algorithm just counts the white pixels.

But in order to get this straightforward algorithm I have to be very good at designing the lighting/optical/mechanical part of the system and maybe that part will cost more than the software.

Maybe in the past I read a sentence like "do as much as possible with the mechanics and as little as possible with the software"; it seems to me it was in a book of the 1990s (or 1980s) about practical machine vision but I can't find the proper citation/reference.

I have found some "proverbs" like:

  Never use software to compensate for a poor lighting system.
  It is not cost effective and will result in a poor system design.

  It is cheaper to add a light-proof shroud to keep sun-light away from
  the object under inspection than to modify the software. Another
  universal truth which is often forgotten.

  Nothing exceeds the speed of light. Any processing that can be done
  optically will save a lot of computer processing later.

in the book "Intelligent Vision Systems for Industry" by Bruce G. Batchelor and Paul F. Whelan and also in B.G. Batchelor and P.F. Whelan (1994), "Machine vision systems: Proverbs, principles, prejudices and priorities", Proceedings of the SPIE - The International Society for Optical Engineering, Vol. 2347 - Machine Vision Applications, Architectures, and Systems Integration III, Boston (USA), pp 374-
383. (see here http://elm.eeng.dcu.ie/~whelanp/proverbs/proverbs.pdf).

The proverbs are also in the 2012 book "Machine Vision Handbook", Editors: Bruce G. Batchelor ISBN: 978-1-84996-168-4.


Otsu binarization with value preservation I'm trying to implement a modified version of the otsu binarization algorithm. I'm trying to binarize document images. But in the binarization procedure I want the object (in this case the text) to retain its original grayscale value while the background takes the value of 255; that is, white. I'm posting a sample image version that I found in a paper.

This is the original image:

This is the resultant image I want to obtain:

Could someone please tell me how to do it in Matlab?

I don't have Matlab handy, but here is how you do it in OpenCV. The example below uses the python interface via Python (x, y):

test = cv2.imread("test.jpg", 0)
(_, otsu) = cv2.threshold(test, 0.0, 255.0, cv2.THRESH_TOZERO_INV + cv2.THRESH_OTSU)
cv2.imshow('otsu', otsu)

This results in your required output:

EDIT : I don't have a copy of Matlab, but I think this is how you would do it (assuming you have the Image Processing Toolbox):

Use graythresh to get the Otsu level, then set anything above that level to white (or 255).

I = imread('doc.jpg');
I = rgb2gray(I);
otsuLevel = graythresh(I);
I(I > otsuLevel) = 255;

Hope that helps!


Is there mathematical relationship between the FFT and DCT transforms I need to model the relationship between the DCT and DFT transforms (If it exists). I mean real signal $x \to y = \textrm{DCT}(x) \to z = \textrm{DFT}(y)$, so I need to get the relationship between the $x$ and $z$ if possible.
For more details:
Assume I have a real signal $x[m],\ m=1,2,….N$, the $N$-point IDCT transform of the signal $x[m]$ is $X[n]$ which can be written as follows (let's ignore the coefficients for simplicity):
$$X[n] = \sum_{m=1}^N x[m] \cos\left( \frac{(2m+1)n\pi}{2N} \right)$$
I need to get the relationship between the FFT of $X[n]\ $ and the signal $x[m]$, so taking the $2N-$point  DFT for the signal $X[n]\ $ (I upsampled the signal to have the DFT of each point resulted from the IDCT), that will give:
$$Y[v] = \sum_{n=1}^{2N} X[n] e^{-\frac{j2\pi mn}{2N}} = \sum_{n=1}^{2N} \left[
\sum_{m=1}^N x[m] \cos\left( \frac{(2m+1)n\pi}{2N} \right) \right] e^{-\frac{j2\pi mv}{2N}}$$
with some mathematical operations, we can get:
$$Y[v] = \sum_{n=1}^{2N} \left( \left[ \sum_{m=1}^N x[m] \cos\left( \frac{(2m+1)n\pi}{2N} \right) \right] \cos\left(\frac{2\pi mv}{2N}\right) - j\left[\sum_{m=1}^N x[m] \cos\left( \frac{(2m+1)n\pi}{2N} \right) \right] \sin\left(\frac{2\pi mv}{2N}\right) \right)$$
I am trying to get the relationship between the $Y[v]$ and $x[m]$ from the relationship above if it exists. How can we get it?

You have some confusion in  the indices, maybe the correct expression would be
$$Y[k] = \sum_{n=1}^N \left[
\sum_{m=1}^N x[m] \cos\left( \frac{(2m+1)n\pi}{2N} \right) \right] e^{-\frac{j2\pi k n}{N}}$$
The DCT is the real part of the odd-indices coefficients of double length DFT
Expressing DCT in terms of complex harmonics
$$z[n] = \sum_{m=1}^N x[m] \left( \exp\left( j \frac{(2m+1)n\pi}{2N} \right) + \exp\left( -j \frac{(2m+1)n\pi}{2N} \right) \right)$$
If you take $X(n)$ being the fourier transform of $x$, in the reals instead of the integers, you could write $z[n] = X(m+1/2)+X(-m-1/2)$, $X(n)$ has Hermitian symmetry and the result agrees with the fact that $z[n]$ is real.
I think that the most efficient way to express the relation between $Y[m]$ and $x[m]$, is that $Y = DFT(DCT(x))$, unless you want to upsample your signal.
The denominator $2N$ in the exponential, the summation is from $1$ up to $N$, and the numerators may be even, so cannot simplify the exponent, the way to go is to compute a $2N$-samples, DFT.
If you create an intermediate signal with odd samples $v_{2i+1} = x_i + x_{N - i}$ and even samples $v_{2i}=0$, and compute $V = DFT^{-1}(v)$, then the first $N$ elements of $V$ equals to $z$, and the $DFT(V) = DFT(DFT^{-1}(v)) = v$


Given a Vertical Derivative Kernel, Find the output image after K is applied to the above input image Im not sure how to solve this, if someone can please help me a bit id appreciate it. I have to solve this by hand. (I know how to solve using matlab).

Read about convolution here. The link has a web tool to understand convolution better. That will help you. Here is a quick primer on how to go about solving your problem.

You have to choose a pixel, place kernel such that centre of kernel is on the pixel you want to compute. Multiply the overlapping pixels and dive by number of non-zero pixels as shown above. For the boundary pixels you have to take a call as to what you will do. There are many options like use zeros as boundary etc. For derivatives it is better to keep the boundaries as zero. Hope this helped you :)


Purpose of image feature detection and matching I'm a new guy in image processing and computer vision, so this question might be stupid to you.

I just learned some feature detection and description algorithms, such as Harris, Hessian, SIFT, SURF, they process images to find out those keypoints and then compute a descriptor for each, the descriptor will be used for feature matching.

I've tried SIFT and SURF, found that they are not so robust as I thought, since for 2 images (one is rotated and affined a little), they don't match the features well, among almost 100 feature points, only 10 matches are good.

So I wonder

What can we use these feature detection algorithms for in practice? Is there any more robust algorithms for feature detection and matching? Or SIFT and SURF is already good, I just need to refine it for further use?
Another problem is that I thought these algorithms are not quite for real-time application (without considering multi-core implementation), but there are some commercial products (such as Kinect) which work and response in real-time! I assume these products also detect and match feature from what they see, do they use the algorithms such as SIFT? How could they detect features so well?
With my limited knowledge, I know feature matching can be used to find out same objects in two images, or estimate homographies, but any other purpose for feature matching?

Image keypoints are a key feature in many Image and Video processing softwares, both industrial and academic. The principle behind is always the same:

detect some meaningful points in some images;
[optional] compute a stable description of the image part surrounding each keypoint;
match keypoints from an image (the template) to another (the query).

Now, some more detail.

Why detecting points? Basically, because a point is an intersection of two lines. As such, it is detected much more accurately (and in a more stable way) as lines or other features (area, etc.). This has been common knowledge in Computer Vision later proved in a paper by Mikolajczyk. This is rooted deeply in any point detector: for example, Harris point detector relies on the Hessian matrix to check that there are two crossing lines.

Why is the description part optional? Basic point detectors (Harris, SUSAN, Moravec, FAST) do not come with a descriptor. The process of matching them is purely location-based. Thus, they are mostly used for video processing. The intuition behind is that you can have many many points, then a robust estimation process (e.g., RANSAC). Recent keypoints however (SIFT and the ones that came after) embed an image patch descriptor, because they are meant to be used in more complex situations (object detection, large baseline matching...).

What's he difference between Harris/Hessian and SIFT/SURF? Besides the descriptor part, SIFT/SURF also include a localization in scale. Thus, when you detectd a SIFT/SURF point, you expect to find it again even if your image has been zoomed in/out, etc.

What's the difference between SIF and SURF? Well, SURF is actually mostly an approximated version of SIFT that is better suited for fast computations (by using integral images). Consequently, SURF is not as stable (as "good") as SIFT under affine transforms, but it comes at something like 1/10th of the computational cost.

Are they suited for realtime software? Yes, fore sure for Harris. SURF was designed to be fast, so a good implementation should be fast too. SIFT was not designed for speed, but you can find GPU implementations. But if you really want to be fast even on lower-end devices, then you need to check recent points such as FAST (detector), BRIEF, ORB, BRISK, FREAK (binary descriptors).

What can you do with them? Well, many things. By tracking points in videos then you can stablize them, augment them with objects, track objects. By matching points between cameras, you can automate the calibration process of 3D reconstruction pipelines. By matching points on parts of object, you can detect objects. And probably much more...


How do window size, sample rate influence FFT pitch estimation? I am trying to create a pitch-detection program which extracts the frequencies of peaks in a power spectrum obtained from an FFT (fftpack). I am extracting the peak frequencies from my spectrum using Quinn's First Estimator to interpolate between bin numbers. This scheme seems to work well under certain conditions. For example, using a rectangular window function with a window size of 1024 and a sample rate of 16000, my algorithm correctly identifies the frequency of a pure A440 tone as 440.06 with a second partial frequency of 880.1. However, under other conditions, it produces inaccurate results. If I change the sample rate (e.g to 8000) or the window size (e.g. to 2048), it still correctly identifies the first partial as 440, but the second partial is somewhere around 892.  The problem becomes even worse for inharmonic tones like those produced by a guitar or piano.

My general question is: In what way do the sample rate, window size, and window function affect frequency estimation of FFT peaks? My assumption was that simply increasing the resolution of the spectrum would increase the accuracy of peak frequency estimation, but this is clearly not my experience (zero padding also does not help).  I am also assuming that the choice of window function will not have much effect because spectral leakage should not change the peak location (though, now that I think about it, spectral leakage could potentially influence the interpolated frequency estimate if the magnitudes of bins adjacent to the peak are artificially increased by leakage from other peaks...).

Any thoughts?

Use a Gaussian window - the Fourier transform of a Gaussian is a Gaussian
Log-scale the spectrum to emphasize peaks and turn the Gaussian peaks into parabolic peaks
Use parabolic interpolation to find the true peaks.

  Note that, as mentioned in §D.1, the Gaussian window transform magnitude is precisely a parabola on a dB scale. As a result, quadratic spectral peak interpolation is exact under the Gaussian window. Of course, we must somehow remove the infinitely long tails of the Gaussian window in practice, but this does not cause much deviation from a parabola, as shown in Fig.3.30.

https://ccrma.stanford.edu/~jos/sasp/Quadratic_Interpolation_Spectral_Peaks.html

I estimate 1000.000004 Hz for a 1000 Hz waveform this way: https://gist.github.com/255291#file_parabolic.py

If you're having trouble, plot the spectrum and use your eyes to see why it's not working.


filtering difference b/w freq and time domain When you want to get the specific frequency (ex. around 1000Hz), you might

take FFT, remain the frequency you want, and reduce the power of the other frequencies you do not.
ex.) if fs=16khz and fftsize=1024, remain the power of bin64 and reduce all other bins
design band pass filtering (Butterworth or others) to pass 1000Hz (with stop/transition band also)

In short, 1 is in frequency domain, and 2 is in time domain.
Then my question is: are these two methods same?
(I think yes in terms of "extracting 1000Hz signal"(Both can do))
If yes, how would you utilize these two methods?
If not, how different are these (and how would you utilize)?
I would appreciate your answers. Thank you.

Multiplication in the frequency domain is the equivalent of convolution (digital filtering) in the time domain - one of the fundamental properties of the Fourier transform.
Using an FFT in practice to filter requires stitching together overlapping chunks and a large delay if in real-time.  However in the frequency domain you can immplement brick-wall filters as easy as any, which is much more costly in time-domain.
The algorithms for fft convolution include "overlap-add" and "overlap-save", these are useful search terms on the matter.
The Python scipy.signal library has an implementation of convolve that automatically selects direct filtering or fft convolution based on an estimate of which will be faster.


The Least Norm Solution of Under Determined Linear System Suppose I have a matrix

$$A=
\begin{pmatrix}
1 &  0  & 1 & 0\\
0  &  1 & 1& 0\\
\end{pmatrix}
$$

where the variables are channel information like assume $X_1$, $X_2$, $X_3$ and $X_4$. Considering the typical least-squares problem

$$Y= AX+b$$

However, matrix $A$ is fat and, thus, $A^{T} A$ is not invertible.  I cannot get $X_3$ or $X_4$ and they will be zero. I can use regularization and optimize joint least-squares and least-norm problem but that does not give me the same $X_1$ and $X_2$. Is there any way to solve this problem?

I tried to stack all the received vectors into matrix $Y$ and do the same for $A$ and $b$ to make matrix $A$ fat. However, due to the structure of $A$, it is still not invertible. Does anyone know how to solve this problem? Can I recover the solution for parameters based on what I get from regularized solution?

Least Squares solution is always well defined for Linear System of Equations.
In your case, which is under determined it means there are many solutions to the Linear Equations.
The Least Squares solution has nice property, it also minimizes the $ {L}_{2} $ norm of the solution (Least Norm Solution) hence it is well defined.
In practice, in order to solve this equation you need to use the pinv() function in MATLAB.
This function uses the Moore Penrose Pseudo Inverse Matrix of $ A $ to solve the equation:
x = pinv(A) * b

It basically decompose the Matrix using the Singular Value Decomposition (SVD) such that $ U S {V}^{T} = A $ and the rebuild the matrix $ {A}^{\dagger} = U {S}^{\dagger} {V}^{T} $ where $ {S}^{\dagger} $ is the result of inverting any non zero element of the diagonal matrix $ S $.
Update
The user @percusse said the following:

You just need A\b
You shouldn't use pinv for linear system solutions that's what I mean. Not only inefficient but also not accurate. There is no unique solution hence you can't talk about a least.
Backslash operator is optimized for all possible linear eq sets. Try it.
There is no least squares solutions for underdetermined systems because of the nontrivial nullspace. Anyways nevermind.

So, to make things clear.
The Least Squares solution always exists and it is also the Least Norm solution.
Namely, from the set $ \mathcal{B} = \left\{ x \mid \min_{x} \left\| A x - b \right\| \right\} $ there is one unique $ \hat{x} $ such that $ \hat{x} = \arg \min_{x} \left\| x \right\|_{2} \text{  subject to  } \left\{ x \mid x \in \mathcal{B} \right\} $.
For instance:
for ii = 1:1e12
    numRows = 3;
    numCols = 27; %<! Must be bigger than 'numRows' for Under Determined System

    mA   = randn([numRows, numCols]);
    vB   = randn([numRows, 1]);
    vXLs = pinv(mA) * vB;
    vX   = mA \ vB;

    if(norm(vXLs) >= norm(vX))
        keyboard();
    end

    if(mod(ii, 1e4) == 0)
        disp(['Iteration Number - ', num2str(ii)]);
    end

end

Remarks:

The difference arises for Under Determined System. For over / well determined system the solutions collide (Namely Tall Matrices which are the "Classic" case in LS problems).
The Psuedo Inverse is much slower to use then the \ operator. If you're looking for any solution, use \. If you're looking for the Least Squares solution (The Least Squares solution is always the least norm solution) for Under Determined System you need to use the pinv() function.
Proof and deeper analysis can be found in my Singular Value Decomposition (SVD) Presentation. Specifically this issue is on pages 50-60.
Cleve Moler's (Inventor of MATLAB) post on Pseudo Inverse vs. \ - The World’s Simplest Impossible Problem.


When is a network called end-to-end training? In machine learning, we often see the expression "end-to-end" learning (or training). However, I do not know that it means. When is a network called end-to-end training? How to recognize a network is end-to-end learning?

From feature extraction to learning the desired result, deep learning algorithms can act as full pipelines for solving tasks at hand. End-to-end learning usually refers to omitting any hand-crafted intermediary algorithms and directly learning the solution of a given problem from the sampled dataset. This could involve concatenation of different networks such as multiple CNNs and LSTMs, which are trained simultaneously.

For the OCR example, instead of trying to classify characters and clustering them into words, it is shown to be a better approach to directly use CNNs to regress the words themselves.

In the optical flow, one could go directly from images to the final flow field, omitting the image difference computations.

In self driving cars, the network can be trained to directly learn how to drive.

In all such examples, the idea is to let the network go from the "raw-est" possible data to the final-most output. This is found to perform better. End-to-end learning reduces the effort of human design and performs better in most applications.


You have got exactly what you requested in your code.

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

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


Implementing a basic demodulator I'm trying to implement a simple demodulator as an evening project (to understand how this works - my background is in physics not engineering) and was wondering if anyone has any advice.

To make a working demodulator I think I need to implement an:

Equaliser (I was thinking of using a steepest descent algorithm).
Channeliser (I don't know how to go about this)
A carrier (or phase-locked?) loop

Does anyone know of any (hopefully free) resources that will help me?

What type of signal do you need to demodulate (BPSK, FSK, OFDM, etc.)? While a lot of practical systems have them, you don't always need an equalizer or phase-locked loop, especially if you're just making a toy simulation where you control all the conditions. If you've never done this before, I wouldn't recommend jumping in with all of that included. Each of those topics that you brought up is a large area of study in itself, so I can't really cover all of the theory in a single answer. I would operate in the following order:

Generate a modulated signal with no noise. Build your receiver and show that you can extract the information data out successfully (since there are no real impairments, you should be able to do this perfectly in simulation).
Add white Gaussian noise to the input signal that your receiver operates on. Test the receiver to calculate its bit error rate at its output as a function of the signal-to-noise ratio in the receiver input signal. Compare these results to theoretical results found in a digital communications textbook (like this or this) to verify that you're operating correctly.
Add a frequency and/or phase offset to the receiver's input signal. Experiment with different techniques for correcting the offset in your receiver (as you noted, a phase-locked loop is one such way). Compare the receiver's bit error rate with and without the phase/frequency synchronizer in place.
Add channel effects to your simulation (i.e. multipath). This can be simulated by convolving the transmitted signal with a simulated channel impulse response. There are a lot of different channel models out there, depending upon the type of scenario that you want to test. At the receiver, add an equalizer that attempts to correct for the channel effects. The easiest-to-implement adaptive equalizer is probably one based on the LMS algorithm. Compare the bit error rate at the receiver with and without the equalizer in place.

This is very high-level and each step is complicated in and of itself, but you're certain to run into trouble if you try to do everything at once with no background.


From Fourier transform to Laplace Transform It's well known that you can estimate the Fourier Transform $X(f)$ of a signal $x(t)$ via its Laplace Transform $X(s)$, just by setting $s = j2\pi f$ to the latter, as long as the region of convergence includes the imaginary axis.

However, I do not have a clear view of how (and when) we can obtain the Laplace Transform via the Fourier Transform of a signal (which is the opposite of what I've stated before).

For example, $x(t) = e^{-at}u(t)$ has a Fourier Transform $X(f) = \frac{1}{a+j2\pi f}$ as long as $a > 0$. Its Laplace Transform is $X(s) = \frac{1}{a+s}$, for any value of $a$, as long as $\mathrm{Re}\{s\} > -a$. We can see that if we set $j2\pi f = s$ to the Fourier Transform, we can directly obtain the Laplace Transform. The same holds for any rational function of $j2\pi f$.

Is there a theorem or something that can be clearly stated about it? It looks to me that it has something to do with the convergence of the Fourier integral (if it does converge, then the Laplace Transform converges as well).

So it sounds like X(s) is the analytic continuation of X(f), given the proper restrictions, is that right?

that's what i would say.  be careful about expressing the function and its argument.
use a different "$X(\cdot)$" if you need to.  like:
$$ X_\mathscr{L}(j 2 \pi f) = X_\mathscr{F}(f) $$
where
$$ X_\mathscr{L}(s) \triangleq \mathscr{L}\Big\{ x(t) \Big\} = \int\limits_{-\infty}^{\infty} x(t) \, e^{-st} \, \mathrm{d}t $$
and
$$ X_\mathscr{F}(f) \triangleq \mathscr{F}\Big\{ x(t) \Big\} = \int\limits_{-\infty}^{\infty} x(t) \, e^{-j2\pi f t} \, \mathrm{d}t $$


How do I find the transfer function in the frequency domain? I was doing some exercises with transfer functions, they were always under the form of $H(z)$ and $H(e^{jw})$ for the frequency response. Today I have found one with $H(f)$. I would like to ask if my solution to the problem is good ?
Let $y(n) = -x(n)$ where $y(n)$ is the output and $x(n)$ is the input. Find the transfer function $H(f)$.

Here is my approach to the problem:

$y(n) = -x(n) => H(n) = \frac{y(n)}{x(n)} = -1$

We apply the Fourier Transform

$H(f) = F\{H\} = F\{-1\}= \int_{-\infty}^{\infty}e^{j2\pi f t}dt = -\delta(f)$

Since $\delta(f) = 1$ we conclude that $H(f) = -1$

Your final answer looks good, but the derivation and notation are a bit mixed.  Remember the capitalized $H()$ refers to the frequency domain transfer function, so writing $H(n)$ with the time-domain index $n$ doesn't make sense.  Also, it is the time-domain impulse response that is a delta function, not the frequency domain transfer function, so $H(f)=-\delta(f)$ is not right.

The problem can be solved very easily:

$y(n)=-x(n)$

Take the Fourier Transform of both sides:

$Y(f)=-X(f)$

Divide both sides by $X(f)$:

$\frac{Y(f)}{X(f)} = H(f) = -1$

In the time domain:

$h(n)=-\delta(n)$


Compression ratio for 4x4 block Say I have the following 4x4 block:

[ 120  121  125  119
  120  120  118  117
  122  121  117  118
  122  122  120  120]

If this block is intra coded the resulting DCT coefficients are the following and rounding to nearest integer:

[ 480.8  3.8  -1.5  0
  -0.1  -1.8  -2.7  2.7
    4   -2.4  -2    2
   1.1   0.2   0    0.3]

Assuming I used QP size of 2 for quantization:

[ 240  2  -1  0
   0  -1  -1  1
   2  -1  -1  1
   1   0   0  0]

How do I calculate the compression ration for the following assuming a pixel is encoded using 8bits/pixel and the DCT coefficient uses 10 bits

I normally calculate the compression by first calculating the

uncompressed bit rate = image width × image height × colour depth × bits per pixels

and then solve uncompressed size/compressed size

But I am not sure how to calculate the compression ratio in this case where DCT coefficient is given!!

First, motivations for an engineer solution, which might be the initial reason of such a question:

The compression ratio of a JPEG encoded image result from a sort of average from different image blocks "compressibility". The JPEG principles rely on two combined aspects:

the average luminosity from block to block varies slowly, hence it can be nicely predicted from preceding blocks. This is the differential prediction used on DC coefficients (here your 240);
local stationarized parts of images can be faithfully approximated by a few quantized 2D cosines.

The scheme you provide is a  simplified version, here are a few. First, the DCT is applied on $8\times8$ blocks, not $4\times4$. Hence, it is more difficult to generate long runs of zeros after quantization. Second, the quantization scheme is linear (dividing by 2). Third, with only one block, you cannot gain from inter-block prediction on DC. Fourth, with  only 15 AC coefficients, the (value,run) Huffman coding might be a little inefficient. To add a few, non-adaptive tables, lack of header, may highly impact the estimated compression ratio  on such a small file.

Last, the "DCT coefficient uses 10 bits" is a bit awkward. So, I interpret it as "The DC coefficient uses 10 bits". With that, you are left with 4 values and their occurrences:

$2$ (2 times)
$1$ (3 times)
$-1$ (5 times)
$0$ (5 times)

from with you can build a Huffman tree, and a non-prefix code-word for each of them. Adding their overall length to the above 10 bits can give you a rough estimate of the total number of bytes actually needed. You can use for instance an  Online Huffman coding tool.

Second, you can work out your own JPEG format, on $4\times4$ blocks, and input an sufficiency large image with those repeated $4\times4$  blocks, to compute actual compression ratios.


Can we break the Shannon capacity? I have a friend working in wireless communications research. He told me that we can transmit more than one symbol in a given slot using one frequency (of course we can decode them at the receiver).

The technique as he said uses a new modulation scheme. Therefore if one transmitting node transmits to one receiving node over a wireless channel and using one antenna at each node, the technique can transmit two symbols at one slot over one frequency.

I am not asking about this technique and I do not know whether it is correct or not but I want to know if one can do this or not? Is this even possible? Can the Shannon limit be broken? Can we prove the impossibility of such technique mathematically?
Other thing I want to know, if this technique is correct what are the consequences? For example what would such technique imply for the famous open problem of the interference channel?

Any suggestions please? Any reference is appreciated.

Most certainly not. While there has been some claims to break Shannon here and there, it usually turned out that the Shannon theorem was just applied in the wrong way. I've yet to see any such claim to actually prove true.

There are some methods known that allow for transmission of multiple data streams at the same time on the same frequency. The MIMO principle employs spatial diversity to achieve that. Comparing a MIMO transmission in a scenario that offers high diversity with the Shannon limit for a SISO transmission in an otherwise similar scenario might actually imply that the MIMO transmission breaks Shannon. Yet, when you write down the Shannon limit correctly for the MIMO transmission, you again see that it still holds.

Another technique to transmit on the same frequency at the same time in the same area would be CDMA (Code Division Multiple Access). Here, the individual signals are multiplied with a set of orthogonal codes so that they can be (perfectly in the ideal case) separated again at the receiver. But multiplying the signal with the orthogonal code will also spread its bandwidth. In the end, each signal employs much more bandwidth than it needs and I've never seen an example where the sum of the rates was higher than Shannon for the whole bandwidth.

While you can never be sure that breaking Shannon is actually impossible, it is a very fundamental law that stood the test of time for a long time. Anyone claiming to break Shannon has most likely made a mistake. There needs to be overwhelming proof for such a claim to be accepted.

On the other hand, transmitting two signals on the same frequency at the same time in the same area is easily possible using the correct method. This is by no means an implication that Shannon is broken.


What are characteristics of shadows in an image given an outdoor setting from an image processing point of view? I'm wondering what basic characteristics of shadows are. I'm looking from a camera image/video processing point of view. I'm considering the case of an outside scene in direct sunlight.

These are the characteristics I could quickly come up with, but I guess there are more.

There are strong edges on the shadow borders
The shadow is ‘darker’ than the surrounding area
The shadow usually ‘deforms’ when an object is rotated
Due to absence of the primary lighting in the shade, the color tends to be dominated mostly by the main diffuse lighting source (which is in most cases the sky --> blue)

I'm wondering which other characteristics I might miss.

The edges of the shadows are crisper where they are nearer to the shadow caster.  Also, since the sun is so far away, beams of light are essentially parallel when they reach us, which means that shadows from sunlight are basically orthographic projections of the shadow casters.  Also, if you knew an objects dimensions and saw its shadow on something that you knew the slope of, you could probably calculate the angle of the sun, and so maybe could calculate what time it was, and maybe something about the season or latitude??


Time frequency analysis for a spreading signal I am a beginner in digital communications: I am studying the spread spectrum communication and I have a question on the spreading signals. For example I have 2 spreading signals and I do a time frequency analysis.

Should the spreading signals overlap in time? And if yes, why?

The overlap is not important, it is just necessary that the spreading sequences  are pseudo-orthogonal.

Can someone confirm?


Does convolution violate duality of a Fourier Transform Duality says that
if
$$
x(t) \iff X(\omega)
$$
then
$$
X(t) \iff 2\pi \cdot x(-\omega)
$$

Given that
$$
f(t)\circledast g(t) \iff F(\omega)G(\omega)
$$
By duality, shouldn't that imply that
$$
f(t)g(t) \iff 2\pi \cdot F(\omega) \circledast G(\omega)
$$
since
$$
f(-t) \circledast g(-t) \cdot 2 \pi = 2 \pi \cdot f(t) \circledast g(t)
$$

When I look at the theorem in my textbook it states that
$$
f(t)g(t) \iff \tfrac{1}{2\pi} \cdot F(\omega) \circledast G(\omega)
$$

I'm confused as to why the constant  is divided by $2\pi$ instead of multiplied by $2\pi$.  Does that mean that duality is violated?

Try using this definition of the continuous Fourier Transform and its inverse:

$$ \mathscr{F} \Big\{ x(t) \Big\} \triangleq X(f) \triangleq \int\limits_{-\infty}^{\infty} x(t) \, e^{-j2 \pi f t} \ \mathrm{d}t $$

$$ \mathscr{F}^{-1} \Big\{ X(f) \Big\} \triangleq x(t) = \int\limits_{-\infty}^{\infty} X(f) \, e^{j2 \pi f t} \ \mathrm{d}f $$

see how well that works for your duality.

Okay, I'll be more explicit.  Like your "$\omega$" definition, we have linearity:

$$ \mathscr{F} \Big\{ x(t) + y(t) \Big\} = \mathscr{F} \Big\{ x(t) \Big\} + \mathscr{F} \Big\{ y(t) \Big\} = X(f) + Y(f) $$

and derivative:

$$ \mathscr{F} \Big\{ \tfrac{\mathrm{d}}{\mathrm{d}t} x(t) \Big\} = j2\pi f \cdot X(f) $$

and delay:

$$ \mathscr{F} \Big\{ x(t-\tau) \Big\} = e^{-j2\pi f \tau} \cdot X(f) $$

But we get duality as:

$$\mathscr{F} \Big\{ X(t) \Big\} = x(-f)$$
$$\mathscr{F}^{-1} \Big\{ x(f) \Big\} = X(-t)$$

We get for convolution:

$$ \mathscr{F} \Big\{ h(t) \circledast x(t) \Big\} = H(f) \cdot X(f) $$

and it's the same in reverse:

$$ \mathscr{F}^{-1} \Big\{ W(f) \circledast X(f) \Big\} = w(t) \cdot x(t) $$

The other cool thing about this unitary definition of the Fourier Transform is how Parseval's theorem come out:

$$ \int\limits_{-\infty}^{\infty} \big| x(t) \big|^2 \ \mathrm{d}t = \int\limits_{-\infty}^{\infty} \big| X(f) \big|^2 \ \mathrm{d}f $$

And the more general Plancherel's theorem:

$$ \int\limits_{-\infty}^{\infty} x(t) \big(y(t)\big)^* \ \mathrm{d}t = \int\limits_{-\infty}^{\infty}  X(f) \big(Y(f)\big)^* \ \mathrm{d}f $$

where $(\cdot)^*$ denotes the complex conjugate.

No funky scaling factors.

This is, in my opinion, clearly the most natural definition of the Fourier Transform if you are never relating it directly to the Laplace Transform.  If you are bopping between the F.T. and the L.T. maybe your definition is better, but then I would say "$X(j \omega)$" so that the F.T. is simply identical to the L.T. but with the substitution of $s=j\omega$.


find the maximum of the DFT of sampled signal The signal $$x(t)=\cos(10\pi t+\phi)+\cos(20\pi t)$$ is sampled with a sampling frequency $F_s$ as $25 \mathrm{Hz}$ where the phase $\phi$ is unknown.

Sampling the continuous time signal $$y[n]=x(nT)$$ where $T=1/F_s$.

After calculating the DFT of $y[n]$ with the length of DFT $N=100$, which of $k$ will give us the maximum of the magnitude of $Y[k]$?

Note that the frequencies of the FFT grid are given by

$$f_k=\frac{kf_s}{N},\qquad k=0,1,.\ldots,N-1\tag{1}$$

where $N$ is the FFT length.

Now observe that the two frequencies of the given signal fall exactly on the grid. So you only have to determine these two frequencies, and then use $(1)$ to figure out the corresponding indices.

EDIT:

It's always good to check your results using some simulation software. The following Octave/Matlab script confirms what you came up with.

phi = 0;    % exact value irrelevant
n = 0:99;
T = 1/25;
y = cos( 10*pi*n*T+phi ) + cos( 20*pi*n*T );
Y = fft(y);
subplot(2,1,1), plot(n,abs(Y),'x'), grid on


Can someone explain complex mapping as appears in Ogata's textbook Refer to diagram above, in Ogata's text on discrete time control, he showed that you can map a curve in the S plane, namely curve 1,2,3,4,5 onto a circle in the Z plane through the complex mapping $e^{Ts}$

What I do not understand is why are we missing the negative branch? This is curious because we must have poles in the z-plane located on the negative branch. So if that branch is missing, how are we suppose to map those poles back to the S plane?

Can someone explain why the negative axis is missing in this mapping?

If we denote $s=\sigma+j\omega$, then for the mapping $z=e^{sT}$, the negative real axis in the $z$-plane in the range $z\in[0,-1]$ corresponds to $\omega=\pm\omega_s/2$ (i.e. Nyquist) and $\sigma\in [-\infty,0]$, where $\omega_s$ is the sampling frequency in radians. So if the curve in the $s$-plane moved along the lines $s=\pm j\omega_s/2$ then the curve in the $z$-plane would move on the negative real axis of the $z$-plane. If it doesn't (as suggested in the figure) then this range of the $z$-plane will be excluded. So if with this mapping you want poles on the negative real axis of the $z$-plane, you need corresponding poles with imaginary parts $\pm \omega_s/2$ in the $s$-plane.


How can I upsample 22 kHz speech audio recording to 44 kHz, maybe using AI? I have a 22 kHz mono audio recording, which is mainly speech, a reading. I would like to upsample somehow to 44 kHz, to improve the audible quality. I have read about there are AI methods for upsampling pictures, even videos to a higher resolution. Maybe there is some similar methods for audio too.
The recording comes from radio broadcast, so the recording was excellent studio quality, streamed in 128 kbps, 44 kHz, stereo, but somebody, who did the recording, downsampled it from the original 44 kHz broadcast to 22 kHz, mono. Hence the sound quality is far from optimal, sounds like an old telephone, I think, due to the missing high frequencies.
The sound sample https://drive.google.com/file/d/1AqdQomHjpKqYM4Wq8SdHGOPDR8Irzjj3/view?usp=sharing

What you want to do is called "Super Resolution" in the field of imaging.
This is an ill posed problem.
Hence in order to solve it you need some prior / model about your audio data.
For example you may look at the model in Speech Super Resolution Using Parallel WaveNet.


time invariance concept? First of all ,please let me know?Is cosine time invariant or time variant?
If latter is the case then if a time varying input like cos is given to a time invariant system,how will it behave?

I have a system $y(t) = x^2(t)$
Apparently from equation it seems to be time invariant
But when i give it input $x(t) = \cos(2 \pi t)$
I do not get time invariant response as shown through following matlab code
Although people are commenting that my system is still time invariant but as shown highlighted in attached photo y(t-3) span is between 3 to 13 while S[x(t-3)] span is between 0 to 10 but as by book both of them should have same time span?I am confused in this difference in time span

clc
clear all
close all
t=0:.001:10;
x=cos(2*pi*t);
y=x.^2;
subplot(3,1,1)
plot(t,y)

legend('y(t)')

subplot(3,1,2)
plot(t+3,y)

legend('y(t-3)')

xn=cos(2*pi*(t-3));
y2=xn.^2;
subplot(3,1,3)
plot(t,y2)

legend('S[x(t-3)]')

You do get a time-invariant response. Your code produces the same output for all three signals. In particular, it produces the same output for $y(\sigma^T\{x(t)\}))$ as for $\sigma^T\{y(x(t))\}$ (plots 2 and 3 in your code). It is hard to see in your case because you have shifted the signal $3\cdot 2\pi$ in time. Whats a cosine shifted by $6\pi$? The same cosine. Everything will be clearer if you plot the signals on top of each other and shift by a more informative value.

Try this code:

clc
clear all
close all
t=0:.001:10;
x=cos(2*pi*t);
y=x.^2;
plot(t,y)
hold on;

plot(t+0.25,y)

xn=cos(2*pi*(t-0.25));
y2=xn.^2;
plot(t,y2)

legend('y(t)','y(t-3)','y[x(t-0.25)]')

hold off

You can see in the attached image, the orange and red lines, corresponding to shifted output and output on shifted input respectively, align with each other. That is why you can not see the red line up until the orange line ends. The unshifted output of the unshifted signal in blue does not align with the two other.


What is Frequency Resolution? Im trying to tackle the following problem while still not having a firm idea on what "frequency resolution" means :

  Suppose we sample a continuous time signal with sampling period Ts = 1/2000, and then use a window of length 1000 on the resulting discrete time signal. If we transform it using a 2000 point DFT what would its frequency resolution be ?

Can anyone help me figure this out ?

Edit:
I've come to realize that my definition below of "Frequency Resolution" is completely wrong (as well as OP's question). Frequency resolution is how similar the window function's magnitude in frequency space is to the Dirac delta function. This is because the product of the window and the signal in the time domain becomes convolution in the frequency domain (and a convolution with the Dirac delta function is a sampling which would give perfect frequency resolution) The fatter the mainlobe (quantified by it's variance), and the higher the sidelobes, the worse the frequency resolution. Additionally, Time resolution can be quantified a the variance of the window function in the time domain.

Frequency Resolution is not Bin Resolution/Width. In the graph below notice that the lobes dont get closer (frequency resolution) even though bin width is decreasing.

Credit: Dan Boschen
Frequency resolution is rather a property of the Fourier transform of the rectangular function (i.e. the sinc function).

We must window functions to work with Fourier transforms (even when working theoretically). As a consequence we are always working with $f(t)w(t)$ rather than the function $f(t)$ itself (here $w(t)$ is a rectangular function). By the Convolution theorem the Fourier transform of a windowed function is always a convolution of $\hat{f}$ with $\hat{w}=$ sinc. Notably when $f$ is sinusoidal, $\hat{f}$ will be a Dirac delta function and the convolution will just be a sampling of a sinc function. Thus we periodically lose frequencies completely when windowing, the periodicity of this loss is the frequency resolution.
Since, on windowed functions, the DTFT is a periodic approximation of the CTFT it also acquires these properties.
The confusion arises because when we dont pad zeros to the DFT (i.e. only sample $f(t)w(t)$ where $w(t)=1$), the bin width is equal to the Frequency resolution.

However, we may also pad zeros (i.e. also sample $f(t)w(t)$ where $w(t)=0$) and this results in the DTF better interpolating the DTFT of $f(t)w(t)$. Confer with the first graph.

To see why the Fourier transform of the rectangular function is the a sinc funcion watch this video and consider the winding of the sinusoidal functions (its quite involved though)

To answer OP's example the bin resolution is $$\frac{F_s}{N} = \frac{2000}{2000}=1$$ where $F_s=2000$Hz is the sampling rate, and $N$ the DFT size.
The frequency resolution is what the bin resolution would be if we just sampled in the window (no zero padding)
$$\frac{F_s}{M} = \frac{1}{T}=2$$ where $M$ is the number of samples in the window, $T$ is the duration of the sample, and $F_s=M / T$.


What is meant by isotropic diffusion? I am working on diffusion of medical images, and I have come across the term isotropic diffusion. Can anyone explain it for me?

As Olli has said, Isotropic means the same in all directions. So the diffusion behaves the same regardless of which direction it is diffusing in.

Heat for example diffuses in an isotropic manner. (Well ... unless you have a strange environment with varying conductivity constants or something ...  but that's another story)

The heat equation is

$\frac{\partial u} {\partial t} = k \nabla ^2 u $

Where $U(x,y)$ is heat at point $(x,y)$ and $t$ is time

$k$ is the thermal diffusion constant.

I think it helps to see visually what this looks like

(From wikipedia)

This is actually just like doing Gaussian Blur on an image.

(Gaussian convolution is a solution to the heat equation)

Gaussian blur before and after

Note that the whole image gets the same blur

Anisotropic diffusion on the other hand may diffuse differently in different directions.

For example, perona and malik suggested doing diffusion based on the Laplacian (that's the 'edginess' $\nabla^2$) of that part of the image.

This can be used for performing blur only in certain parts of an image.

See here for an image which has had anisotropic diffusion. This has had the effect of removing noise, while preserving the main features of the image.

The equation proposed by Perona and Mailk is

$\frac{\partial U}{\partial t} = c(x,y,t) \nabla ^2 U + \nabla c . \nabla U$

where $c $ is a function such as

$C(||\nabla U||) = e^{-(||\nabla U||/K)^2}$ for example


Tracking movement of audio speaker driver using images I have a speaker emitting low-frequency sinusoidal wave and I have physically marked a curved edge on the speaker driver as shown in red in the image below, which I would like to track:

I detect the movement of the speaker driver using a Dynamic vision sensor(DVS) based event-based camera which can capture changes in intensity at every pixel with very high resolution(in the order of microseconds).

Can anyone point me towards how to model and predict the movement of this red edge on the speaker driver, moving back and forth in the image? If I were to plot the column displacement (with respect to the image) of this edge versus time, I am guessing the plot would resemble a sinusoid matching the frequency of the sound wave the speaker is emitting. Based on this, I think I can use a Kalman filter to predict the sinusoidal movement but I am not sure how to set up the equations.

I would appreciate any tips, pointers or resources on this.
Thank you!

EDIT:

The speaker is emitting a sinusoidal wave of single, constant frequency and amplitude.

This is the setup of the entire experiment:

A single frequency coming out of speaker is very much a stationary process, then why go for a Kalman filter? Use a Weiner filter instead, much more easier to setup and makes sense for your problem.

Just go for a standard weiner filter, you could go for recursive least squares as well to model this.


Effect of Gaussian Blur on Different Frequency Components of an Image Gaussian blur is lowpass filter which means it pass the low frequency components  and stop the high frequency component but the point  I am not getting is that if it is lowpass filter then how it is possible to prevent the edges which is nothing but the "HIGH FREQUENCY COMPONENT"?

Gaussian Blur is Spatially Invariant Linear Filter.
Hence it can be analyzed in the Frequency Domain which in fact shows its Low Pass properties.
Namely it attenuates High Frequency Energy.

In Image, Edges, which are abrupt change from one color to other, requires high frequencies in order to be local.
Gaussian Blur attenuates and smear it.
As written above since it is Spatially Invariant, it behaves the same on DC zone or Edge.

Now, there are families of "Edge Aware Filters" which are adaptive to the neighborhood they operate at:

Bilateral Filter.
Non Local Means Filter.
Median Filter.
Guided Filter.

All of those react to the properties of the pixels they operate at hence can not be defined in Linear Convolution (Spatially Invariant).


How to measure the frequency response of an acoustic filter by input and output signals I am trying to measure the frequency response of acoustic filter. I sample the sound before propagation in the filter and after by 2 different mics. And now I am trying to find the system (the acoustic path) frequency response.

Using matlab how can I do it?

I tried the identification app but the result is not so reliable and fails if the signal is changed.
I tried $\displaystyle\rm H=\frac{\lvert FFT(output)\rvert}{\lvert FFT(input)\rvert}$ but the result is bad too.

UPDATE: After running this code:

[y, Ss] = audioread('4.wav');

h = fir1(30,0.2,rectwin(31));
y_f = filter(h,1,y);

[Txy,f] = tfestimate(y,y_f);
figure;
plot(f,mag2db(abs(Txy)))

I am getting:

What I doing wrong?

tfestimate is a good approach. But it seems you're using the wrong kind of input to test the filter. If your input signal only contains lower frequencies tfestimate has not enough information to compute the transfer for higher frequencies. It's sort of a division by zero situation, then.


CDMA scheme that preserves constant envelope (NRZ/ASK shape) (This is a followup question to Is this a "valid" CDMA modulation?)

Is it possible to find a CDMA scheme in which the receiver can only have a simple slicer?
In other words, for the receiver, the signal should look like an ASK-like signal, i.e. it can detect a bit of 1 (above the average) or a bit of 0 (below the average) but cannot detect signal levels. It also cannot perform despreading in analog domain.

Further explanations:

It seems to me that any CDMA scheme that uses Walsh codes (or PN codes) result inevitably in a signal with multiple levels when the transmissions of multiple stations are summed. For example, even in the simple example from https://en.wikipedia.org/wiki/Code-division_multiple_access#Example, the output signal is $(0, −2, −2, 0, 2, 0, 2, 0)$, consisting of three different signal levels. For demodulation, the receiver needs to sample the signal levels in order to sum them together.

On the contrary, operating in $GF(2)$ seems to preserve the desired signal properties. However, it seems real spectrum spreading is not possible with this method, since only one station would transmit a bit at a given time (and while technically being CDMA, being more time division multiplex, see Is this a "valid" CDMA modulation?)

Yes you are correct that the reception of multiple CDMA signals from different transmitters tends toward a complex Gaussian distribution given the central limit theorem at play - which is not very different from the thermal noise that a spread spectrum signal may be buried in at reception. Yet this does not preclude us from making a hard decision on each sample at the receiver! So there is no need for a constant envelope CDMA in order to use a simple slicer.

Since CDMA signals are buried in the "noise" of other signals and natural noise, if you take multiple samples and average you can get higher effective precision depending on the number of times you over sample within any symbol interval. Oversampling and averaging in the presence of noise is called "dithering"; as long as there is zero-mean noise present the average will converge to higher precision. (Simple example: consider a 1 bit converter that puts out 0 or 1 based on the input being higher or lower than 0.5. If truth was .25, and if we had a uniform noise distribution on each sample that exceeded +/-0.5, then over many samples we will decide 3 zeros for every 1 zero and the average with be 0.25.)

If you do this simplistically by just increasing your sampling rate, you end up getting an extra $0.5$ bits of precision for every doubling of the sampling rate. See this post for further details: What are advantages of having higher sampling rate of a signal?

Sigma Delta converters do this with noise shaping for even higher gains in precision vs the sampling rate.

Given the reception of CDMA is to correlate by accumulating the results of each sample over long codes before deciding what symbol was actually sent (see How CDMA receiver extract it's corresponding data from the receiving modulated & superposition-ed signal?) , it is not even necessary to oversample to achieve the same effect of averaging by dithering as described above- the noise gets averaged out in the correlation.

This is demonstrated below as applicable to both a standard bi-phase (BPSK) communication and BPSK with CDMA with the descriptions above. In the plot we see a barker code sequence [1 0 1 1 0 1 1 1 0 0 0] which has the desirable property as a PN sequence to not be correlated to shifted versions of itself, which was detailed at the link given above. As transmitted it is either sent as a $+1$ or a $-1$. As received it is corrupted by noise and the histogram of all received signals given a $+1$ was transmitted (blue) or a $-1$ was transmitted (red) is shown on the left. For the samples in the histogram that go beyond the zero-crossing, an error would result when we make a decision on that sample. In standard BPSK communication we are indeed concerned about the errors in every bit, to the amount of that histogram crossing the zero-crossing is a big deal and results in the waterfall curves we get relating BER to SNR (as noise in the SNR is the standard deviation of that histogram related to S which is the average).

HOWEVER notably in CDMA any individual error is of less concern since we are no longer concerned with the individual decisions (hence we will call the individual values in the sequence chips and not bits) but the overall average of our decisions over the complete symbol interval.

Consider the two plots below, which is the same me as our plot above showing a BPSK sequence in high and low SNR conditions. The upper plot with what would result at a higher SNR and a relatively lower error rate on each chip. The bottom plot shows the typical CDMA application where we we receive the signal in the presence of very high noise, including the signals from other CDMA transmitters using other codes (so to us they appear very much as AWGN noise). In this case, the tail showing the distribution of samples that would result in an error on each chip is quite large. However to be clear in CDMA we are capturing lots of samples before making a final decision as to which overall code was transmitted--an individual chip error is not a significant consequence. So in the lower plot as long as we are averaging a very large number of samples before we make our final decision (In GPS there are 1023 chips in each sequence, and we also go through 20 sequences before making a data bit decision, so 20460 samples if we only used one sample per chip) we will still demodulate with a very low error rate.

If that is still confusing consider that we undo the code sequence by multiplying it with the code in the receiver, so if a data symbol of one is desired, the transmitter will send the sequence as [1 -1 1 1 -1 1 1 1 -1 -1 -1], and if a data symbol of of zero is desired the transmitter with send the inverse of the code as [-1 1 -1 -1 1 -1 -1 -1 1 1 1]. Thus we convert the code sequence to all one's for data symbol "one" and all minus one's for data symbol "zero". In the presence of AWGN noise only a simple repetition code would have the same effect.

However to note: A typical bi-phase CDMA signal (such as GPS) given the typical phase, time and frequency offset of the receiver will be equally distributed in AM and PM components. Slicing is equivalent to hard-limiting which removes all AM. You may conclude from this that there would be a 3 dB loss of signal; actual results range from 1.6 to 2.2 dB depending on the excess bandwidth of the receiver prior to slicing (more bandwidth less noise penalty). There is a nice chart showing this in James Spilker's GPS books "Global Position Systems Theory and Practice Volume I and II".

This is illustrated further with the graphic below depicting operation of CDMA in the presence of a non-coherent single-tone jammer where in fact there would be a 3 dB loss in SNR under this condition due to hard-limiting the signal. This helps illustrate the important distinction on the effects of hard limiting a CDMA signal that we would expect to receive with a negative SNR, compared to a traditional BPSK signal that we would expect to receive with a positive SNR.

For high SNR conditions with the signal much larger than the jammer, the diagram on the left shows the signal as received with some arbitrary phase offset (and would also have a carrier offset, so this graphic would be rotating at that offset rate). The fixed phasor represents our signal relative to the jammer that is shown in the smaller phasor rotating with its relative offset in frequency (non-coherent). Here the noise due to the jammer is composed of equal AM and PM components, and with hard limiting we remove all the AM without changing our signal level resulting in a +3 dB improvement in SNR.

For low SNR conditions, the situation is reversed in that our signal is much smaller than our noise which is the typcial case for CDMA operation. So here our signal prior to hard-limiting is the AM and PM components on the single tone jammer. With hard limiting we reduce our signal level 3 dB further relative to the jammer, thus is a loss of 3 dB in SNR.

The other issue with hard-limiting and slicing is vulnerability to coherent jammers. If your receiver is slicing and an intentional jammer is able to send a single tone coherently in phase then then received signal will be an AM only component on top of the jamming signal and will be completely removed.  Jamming and interference rejection in general is a high motivation to use higher levels of quantization for CDMA signals; there is little benefit of going beyond 3-4 bits otherwise. So for AJ considerations it is a really bad idea to hard-limit or equivalently slice.