linear convolution toeplitz matrix vs circular convolution toeplitz matrix I have an issue in understanding the difference between building the Toeplitz matrix when the convolution is linear and when it's circular. As I know that Toeplitz matrix $H$ can be built as following
H = toeplitz(h;zeros(N-L,1),h(1), zeros(1,N-1));

where h is the channel and L is the length of the channel, and N is the total length of symbol convoluted with the channel.
My question, if the convolution with the channel is circular, will the Toeplitz matrix still be built in the same way?

They are in general different. For two signals of lengths $N$ and $M$, linear and circular convolution are equivalent if the output is specified to be of length $N + M - 1$ with the appropriate padding. Convolution via the DFT is inherently circular, which is why padding must be done before the inverse DFT to yield the linear convolution. So, this is a special case where they are the same.
If your goal is to always yield linear convolution, then don't worry about forming a circular Toeplitz matrix since the result will be the same when using the regular Toeplitz and is simpler to do so.
Below is some sample code and output where we form regular and circular Toeplitz matrices with a specified output of length $N + M - 1$:
%% Toeplitz Convolution

x = [1 8 3 2 5];
h = [3 4 1];

% Form the row and column vectors for the Toeplitz matrix
r = [h zeros(1, length(x) - 1)];
c = [h(1) zeros(1, length(x) - 1)];

% Toeplitz matrix
hConv = toeplitz(c,r)

% Compare the two types of convolutions
y1 = x*hConv
y2 = conv(x, h)

hConv =

 3     4     1     0     0     0     0
 0     3     4     1     0     0     0
 0     0     3     4     1     0     0
 0     0     0     3     4     1     0
 0     0     0     0     3     4     1

y1 =

 3    28    42    26    26    22     5

y2 =

 3    28    42    26    26    22     5

%% Toeplitz Circular Convolution

% Convolution length
n = length(x) + length(h) - 1;
numElementDiff = n - length(h);

% Set up the circular Toeplitz matrix
c = [h(1) fliplr([h(2:end) zeros(1, numElementDiff)])];
hConvCirc = toeplitz(c, [h zeros(1, numElementDiff)])

% Compare the two types of convolutions
y1 = [x zeros(1, length(c) - length(x))]*hConvCirc
y2 = cconv(x, h, n)

hConvCirc =

     3     4     1     0     0     0     0
     0     3     4     1     0     0     0
     0     0     3     4     1     0     0
     0     0     0     3     4     1     0
     0     0     0     0     3     4     1
     1     0     0     0     0     3     4
     4     1     0     0     0     0     3

y1 =

     3    28    42    26    26    22     5

y2 =

    3.0000   28.0000   42.0000   26.0000   26.0000   22.0000    5.0000

Here we're testing three things:

Linear convolution conv() is equivalent to performing the matrix multiplication with the appropriate Toeplitz matrix.
Circular convolution cconv() is equivalent to performing the matrix multiplication with the appropriate circular Toeplitz matrix.
Output length is specified as  $N + M - 1$, so we see that linear and circular convolution are equivalent.

If you are going to perform circular convolution of varying sizes, then you must form the Toeplitz matrix differently. This usually involves some type of padding with the matrix entries themselves or the signal(s) being operated on. Mathworks has a good summary of Toeplitz matrices here and linear vs circular convolution here.


Is it possible to transmit and/or receive multiple types of signals at the same time? Using an inexpensive SDR (such as the USRP or Per Vices Phi), I want to process two or more communication protocols at the same time (e.g. WiFi and GSM). Is it possible to Tx (or Rx) both at the same time, like FDMA, or would they need to do some fast switching between modulations and frequencies, like TDMA? I'm not sure if the receiver would be able to accept it though since it is not expecting TDMA.

So I think FDMA would be the best approach, but like I said the question is whether or not it's possible to construct a waveform that can transmit both protocols simultaneously. I think it may also depend on frequency separation, combined bandwidths, mutual sample rate, etc.

Theoretically it's possible to do a frequency multiplex of two different protocols with one SDR transmitter. However, the practical main limitation is the speed of the digital-to-analog converter (DAC). E.g. for a WLAN and GSM band the channels that are closest on the frequency axis are at 1800 MHz and 2412 MHz, respectively. For generating such a signal you'd need a sampling rate of over 1200 MSa/s. Even if you find an SDR device that is equipped with an DAC that fast, it will be very challenging to filter the aliasing spectra generated by the DAC in order to avoid out-of-band radiation. So your sampling rate has to be even higher in practice.

If you want to transmit two independent data streams at bands very close to another, OFDM can be advantageous because of it's narrow-band characteristic. If the modulation is synchronized, the bands can even overlap. The latter is an case without practical relevance, though.


How to interpret results of Fourier transform by using Python Connected with this topic , I would like to have some clarifications about the Fourier Transform and the Python tools that make it.
Having my data set whose plot is showed in the above link, my array of data is stored in rate. Then, if I check the result from the Fourier transform:

yfft = np.fft.rfft(rate)
y_smooth = np.fft.irfft(yfft)
print y_smooth

I get:

[  3.20000028e+01+0.j          -7.92633533e+00-7.75032655j
-5.34245688e+00-3.82390808j  -2.39880935e+00-0.60037314j
-8.71249166e-01-0.04362804j  -8.52507194e-02+1.09973189j
-1.58332047e-01+1.15256295j  -3.38909974e-01+0.96327143j
 5.51323015e-02+0.43936749j   3.81205886e-01-0.02489853j
 2.99110153e-01-0.26319101j   3.28283652e-01+0.11668814j
 1.08849205e-02-0.11317937j   3.79952011e-01+0.22278952j
 7.37720998e-02-0.26810083j  -2.14247586e-01-0.25950865j
-1.32242152e-01+0.j        ]

So, I am quite confused about this result, and I would need some tips how to interpret these numbers.

1) The array contains $17$ complex number : does it mean that my signal can be approximated by a superimposition of $17$ harmonics?

2) Are these numbers the power of each harmonic, or what?

3) are they ordered by the frequency order, or by power? Let's say I want to put to zero the higher orders, which ones should I put to zero?

The array contains 17 complex number : does it mean that my signal can
  be approximated by a superimposition of 17 harmonics?

16.

The first component is the DC component.

Note that this is not specific to your signal. Any real sequence of length 32 can be decomposed as 1 real DC component, 16 complex amplitudes, and 1 real amplitude (32 real number ins, 32 real numbers out).

  Are these numbers the power of each harmonic, or what?

They are the complex amplitudes.

  are they ordered by the frequency order, or by power? Let's say I want
  to put to zero the higher orders

They are ordered by frequency.

  which ones should I put to zero?

The last ones, but why would you want to do that?


Transfer function estimation from logarithmic sweep frequency response data I have complex frequency response data (of an analog system) in the range of 100 Hz to 100 GHz, and it is sampled in frequency with logarithmic spacing. I would like to be able to turn this into a filter in MATLAB such that I can multiply it with the fft of time-domain signal.

I'm not seeing a good match when using invfreqs.  What is the best way to do this?

Use FDLS (Frequency Domain Least Squares) to create a model of your measured frequency response. Use the frequency response of that model, evaluated at the FFT bin frequencies, as your filter.

Alternate: if a low-order IIR (or perhaps even FIR) model provides a good-enough fit, then use the model as your filter for direct time-domain convolution.


In opecv function Hough circles how does parameter 1 and 2 affect circle detection and how can I adjust them to increase accuracy? I'm using opencv houghcircle function to detect circles.

http://docs.opencv.org/doc/tutorials/imgproc/imgtrans/hough_circle/hough_circle.html

Problem is if I set minimum radius to null it detects wrong center and draws smaller radius circle or sometimes larger radius circle.  I'm looking for increasing accuracy of detection by adjusting threshold parameter 1 and 2.

Intuition for parameters of HoughCircles:

image: 8-bit, single channel image. If working with a color image, convert to grayscale first.

method: Defines the method to detect circles in images.  Currently, the only implemented method is cv2.HOUGH_GRADIENT, which corresponds to the Yuen et al. paper.

dp: Resolution of the accumulator array.  Votes cast are binned into squares set by dp size.  Set too small and only perfect circles are found, set too high and noise collaborates to vote for non-circles.

minDist: Minimum distance between the center (x, y) coordinates of detected circles. If the minDist is too small, multiple circles in the same neighborhood as the original may be (falsely) detected. If the minDist is too large, then some circles may not be detected at all.

param1: Is a number forwarded to the Canny edge detector (applied to a grayscale image) that represents the threshold1 passed to Canny(...). Canny uses that number to guide edge detection: http://docs.opencv.org/2.4/modules/imgproc/doc/feature_detection.html?highlight=canny

param2: Accumulator threshold value for the cv2.HOUGH_GRADIENT method.  The smaller the threshold is, the more circles will be detected (including false circles). The larger the threshold is, the more circles will potentially be returned.

minRadius: Minimum size of the radius in pixels. Don't set minRadius and maxRadius far apart unless you want all possible circles that might be found in that range.

maxRadius: Maximum size of the radius (in pixels).  Don't set minRadius and maxRadius far apart unless you want all possible circles found in that range.

If you don't know radius, you must auto-tune hyper parameters:

If you don't feed the right radius (+- 15 pixels) to HoughCircles, the output is next to useless because it returns  circles everywhere and nowhere without rhyme or reason.

However there is a trick if you don't know the radius and don't know the optimal dp resolution and accumulator array threshold.  HoughCircles is quick (linear time?), put HoughCircles in a loop, start out with a big minDist, fine dp resolution, and high accumulator voting threshold, and start guessing the largest radius of circle you're willing to find, HoughCircles returns nothing.  Each iteration, slowly decrease radius, decrease threshold, increase dp resolution.  Eventually HoughCircles will find a pixel-perfect circle.

Do this process a second time and approach from the other side, loop HoughCircles starting out guessing a minimum radius, low resolution, high voting threshold, eventually the smallest and highest quality circle will be found.  If the two circles agree on a common circle, you can be very confidant the circle found is high quality, otherwise increment/decrement the hyperparameters at different ranges and speeds.

HoughCircles is mathematically precise and exactly correct, it has been proven to find pixel perfect circles even when you don't know radius.  With hyper-parameters auto-tuned correctly for your images, the circles found are pixel perfect.

Not only can you use HoughCircles to find the middle of the circle, you can use it to define the outer lip and inner lip of the circle's edge gradient (sub pixel accuracy). For example, the outer lip or inner lip of a coin's circumference.

Many passes of many iterations must be used to auto-tune hyper parameters uniquely for every image to slowly discover which circles are the one you want.

Prove pixel perfect accuracy:

If performance must be proved as sub-pixel perfect for an industrial setting, you can't get around intuitively understanding what the algorithm is doing so you can fiddle with its implementation to suit your needs.  For in depth guide, See: "Computer Vision: A Modern Approach" by David Forsyth and Jean Ponce, 2nd Edition, Chapter 10 on the Hough Transform.


What is the difference between convolution filter, low-pass filter, and median filter? Isn't blur filters, like median filter, a type of convolution filter?
How does a low-pass filter relate to them?
Why does this guy differentiate convolution filters here from blur filters here?

A median filter is most certainly not a "blur" filter, purely on the basis that it tends to preserve edges. Edges are abrupt transitions of brightness and therefore that information is encoded in the high frequencies of the spectrum. Incidentally those high frequencies are the ones that low-pass filters suppress the most, leading to that "blurry" appearance because only the form of the depicted objects is retained rather than their details.

Convolution can only be used to represent linear time invariant systems. That is, systems whose relationship between the input and the output is a linear combination and therefore proportional inputs produce proportional outputs. Furthermore, this relationship does not change with respect to time.

A median operator works by sorting the values of the $M \times N$ mask of pixels surrounding some $i,j$ pixel and then assigning their median (the pixel value that happens to lie at the midpoint of the range of values) to the $i,j$ pixel.

This operation, of sorting and then selecting the median, cannot be expressed as a linear combination of the pixel values within the $M \times N$ mask, hence the name "non-linear filter".

Hope this helps.


JPEG: In image compression why are higher frequency components divided by a larger number during quantization Once 2D FDCT has been carried out on 8x8 blocks of the image data, the coefficients are quantized. This involves use of a quantization matrix where we simply divide each FDCT coefficient by the corresponding value in the 8x8 quantization matrix.

The quantization matrix has bigger numbers on the bottom right and smaller numbers on the top left. Bottom right are the highest frequencies and conversely, top left is DC. Why is quantization not uniform?

Note that the matrix which divides the DCT coefficients, performs quantization and the larger the divisor the coarser will be the step size at that frequency.

And that's the very essence of the perceptive image coding, where less perceived (roughly less important) content is coarsely quantized while the more important and dominant low frequency content is finely quantized and represented better. This will result in more quantized zeros at high frequencies which eventually yield more compression while maintaining less distortion compared to a uniform quantization.

Hence this takes advantage of human perception characteristics. Had you thrown away the low frequency information most people would object immediately to the result, while they would not even recognize the loss of high frequency information until some aggressive threshold.


How to reprocess SDR data to be equal to that captured with a higher RBW I have a data set that consists of max hold frequency vs intensity data from an SDR at a specific RBW. I want to reprocess the data as if I sampled it with the larger resolution bandwidth so I can compare what that file might have looked like.
What would be the correct way to go about doing this to get as close to what my SDR's spectrum analyzer would have looked like with the larger RBW.
My approach was going to be using a window function using the size of the larger RBW and recalculating the spectral density in that window.
I am using scipy to process the data.

Max-hold would not be the best data to work with since it represents the peak envelope instead of the average power at each bin. If working with post-detection data captured from a spectrum analyzer, even with max-hold off there are additional but smaller correction factors which are only of concern if an accurate noise measurement is desired. For those cases, I linked an app-note with further details at the very bottom of this post.
Ideally the OP has access to the time domain data, in which case windowing the data would result in the equivalent of a wider resolution bandwidth (RBW), and be a simple approach for proceeding with reasonable ranges for the increase in RBW. If not, the equivalent operation on the frequency samples can be done using a circular convolution with the Kernel of the time domain window (the Discrete-Time Fourier Transform of the window).
A Gaussian window would likely approximate that used for most spectrum analyzer measurements, but a DPSS window would be optimum for minimizing time-bandwidth product (as desired for spectrum analyzers) given finite measurement intervals, or the Kaiser window which is very close to the DPSS minimum. For further details on that see this post on emulating a Spectrum Analyzer and detailing the window choices further:
Classical spectrum analyzer model
These other related posts on resolution bandwidth and equivalent noise bandwidth with windowing may be of additional interest:
Find the Equivalent Noise Bandwidth
How to calculate resolution of DFT with Hamming/Hann window?
Additionally the SDR's assumed resolution bandwidth and the increase in RBW can be confirmed by computing the autocorrelation using the frequency results before and after the additional processing, assuming the captured signal itself is reasonably white noise within the sampling bandwidth.
See page 6-8 of this app note (App Note 1303, old HP)  for further details on additional corrections up to 2.5 dB when using the log power or lin power measurements from a spectrum analyzer in order to get "true-rms" power measurements assuming noise floors as limited by Additive White Gaussian Noise (AWGN).


Autocorrelation of the product of deterministic and random signal I was wondering how to calculate the autocorrelation of a deterministic signal $x(t)$ multiplied by a stochastic process $M(t)$, whose autocorellation $R_M(\tau)$ is known a priori. In my case, $x(t)$ is a truncated monolateral exponentially decaying function.

I suppose that the result of such multiplication  $y(t) = x(t) \cdot M(t)$ is again a stochastic process, but when approaching the calculation of the autocorrelation of $y(t)$ I obtain something that is not even, therefore I suppose I am making some mistakes. I know that the definition of autocorrelation for deterministic signals is different from the one of stochastic processes, but I do not know how to connect the two of them.

As correctly pointed out in the comments, in general the process $Y(t)=x(t)M(t)$ is not wide-sense stationary (WSS), i.e. its autocorrelation function depends not only on the time difference parameter $\tau$, but also on the absolute time $t$:

$$R_Y(\tau,t)=E[Y(t+\tau)Y^*(t)]=E[x(t+\tau)M(t+\tau)x^*(t)M^*(t)]=\\
=x(t+\tau)x^*(t)E[M(t+\tau)M^*(t)]=x(t+\tau)x^*(t)R_M(\tau)\tag{1}$$

In your case you just need to evaluate (1) with the given function $x(t)$. However, as expected you will end up with a function of two variables because $Y(t)$ is not WSS.

There are a few special cases in which the resulting process is indeed WSS or can be easily made WSS. The first case is modulation by a complex exponential:

$$x(t)=e^{j\omega_ct}$$

in which case

$$x(t+\tau)x^*(t)=e^{j\omega_c(t+\tau)}e^{-j\omega_ct}=e^{j\omega_c\tau}$$

only depends on $\tau$ and not on $t$. Another case of interest is the case where $x(t)$ is $T$-periodic. In this case the process $x(t)M(t)$ can be made WSS by introducing a random phase epoch $\Theta$ with uniform distribution in the interval $[0,T]$ which is independent of $M(t)$:

$$Y(t)=x(t+\Theta)M(t)$$

The autocorrelation function of $Y(t)$ is then given by

$$R_Y(\tau)=R_M(\tau)\frac{1}{T}\int_0^Tx(\alpha+\tau)x^*(\alpha)d\alpha$$


Why LFO should be alias-free? I'm really tripping about this point! I try to reason with you and explain what's my trouble.
Consider a scenario with Sample Rate 44100hz and a LFO square wave that will modulate some internal parameter of my plugin.

Its clear that if I generate that signal using this function:

if (mPhase <= mPI) {
    value = 1.0;
}
else {
    value = -1.0;
}

mPhase += mPhaseIncrement;
while (mPhase >= mTwoPI) {
   mPhase -= mTwoPI;
}

once it reaches my D/A converter, the interpolation would make a signal that will contains alias (the more I pitch up foundamental, the more alias I got).

That's normal: I don't have sufficient samples to represent all partials of the signal I'd like to represent. Clear.

But, why this is true also on LFO?

I mean: once I create my digital square (that, until it interpolate, it's really just 1/-1), the target param would be modulated really as a "real" square wave.

Why it should modulate the knob in a way that contains "alias"? The knob won't move as the signal being generated after the interpolation (with ripples and not linear movement). It really moves linear as a square wave: -1, 1, -1, 1. No "alias".

I hope you got what I mean? I can't see visually the "alias" introduced by LFO (by square wave, in this case) modulating a knob.

Note: I know that using an LFO with the code above at a frequency lower than 20hz produce "un-noticeable" alias for audio/human hearing. But for the reason I question, it shouldn't produce alias at all.

TLDR; The LFO should not be anti-aliased. The final signal, after applying the modulation, should be alias free. In practice, you either ignore this or approximate the latter.

Details:

The correct (although not necessarily always practical) way to think about this is to think of the whole signal processing/synthesis chain as a continuous time process, and then apply an anti-aliasing filter in the end, before discretizing.

So say, for example, that you're modulating the amplitude of an audio oscillator with a square LFO, i.e. $x(t) = (1 + 0.5 v_\mathrm{LFO}(t))v_\mathrm{osc}(t)$, where $v_\mathrm{LFO}$ is the output of the LFO with amplitude $1$ and $v_\mathrm{osc}$ is the output of the oscillator, and assume that the audio oscillator is perfectly anti-aliased. Now, the x(t) will have discontinuous edges at transitions of the LFO, which will generate an infinite number of frequency components, which will alias. However, if you had done the modulation in the continuous time domain and then sampled with appropriate anti-aliasing filters, the transition edges would be smoothed and the signal as a whole would be band-limited.

Another effect is that the spectral components of the LFO itself will be aliased, which may or may not have audible significance. This can be seen for example by the fact that if the period of the LFO is not an integer number of samples, the lengths of the sequences of $1$'s and $-1$'s will vary (by one sample) depending on which side of a sample instant the transition happens to land on on which cycle.

As a further extreme example, say you're modulating the amplitude of a constant DC offset. In this case the output will be that of the LFO itself, and if you bring it to audio frequencies you will certainly hear aliasing if the LFO isn't anti-aliased.

On the other hand, it is not correct to apply an anti-aliased LFO to an anti-aliased oscillator, since modulation is a non-linear effect: ideally you should apply a continuous time LFO to a continuous time oscillator, and anti-alias (i.e. sample with bandlimiting) the result. However, smoothing the LFO transitions a little may restrict the bandwidth of the frequencies produced by the modulation, thus partially remedying the extra aliasing. You should not use a sharp anti-alias filter though, since its step response ripple could cause unwanted modulation, but rather a filter which has a good phase response (simple one-pole or a higher order Bessel).

As far as I know, there is no general method to remove aliasing in a system with modulation, and if there were, it would certainly depend on the specific modulation and what is being modulated (frequency vs. amplitude, oscillator vs. filter, and so on). However, oversampling is a straightforward if computationally costly method which can be used to reduce all aliasing in the audio band, including that caused by modulation.

Finally, as you say, if the LFO is well below audio frequency, it is quite unlikely that any of this will be audible. However, you should be careful in that also modulation depth will have and effect: if you're for example modulating the frequency of an oscillator, even with a low frequency, but such a large modulation depth that the change in the frequency is appreciable over a period of the oscillator, you may get audible effects (at least in principle, I have not made experiments nor do I have a reference to cite).

Addendum:

Why do we care about the continuous time model on a discrete system, when working with audio?

(I'm of course assuming here you mean by plug-in a VST or some other audio processing plugin).

Consider what is your end goal here: to produce a sound to be heard by a person. That is inevitably a continuous time (CT) process. What you have control over is the stream of samples sent to the DAC, which is where the signal is transformed from discrete to continuous time.

In other words, you are actually in the business of producing continuous time signals, and the discrete time (DT) representation is simply a tool to reach that end goal. This is what I mean by talking about discontinuities: they are discontinuities in the imaginary CT signal, samples of which form your DT signal.

The easiest way then to reason about the output CT signal implied by your DT signal is to form the DT signal by properly sampling a mathematically formed CT signal. For example with a sawtooth oscillator, you take an ideal mathematical ramp, you bandlimit it and sample it (by whatever method you like, such as BLIP, BLEP, oversampling etc), and then the sampling theorem will tell you that the output of the DAC will be that same bandlimited sawtooth. If you had instead thought of the sawtooth as a DT process and simply formed a naive discrete ramp, you would have been gravely disappointed, as you would have gotten the nasty sound of an aliased sawtooth, because the DT process you used happens to be equivalent to the naive sampling of a CT non-bandlimited sawtooth. I.e. your ears agree with the conclusions you get from thinking of the DT signal as a sampling of a CT signal!

Note that the above applies just as well to all parts of the synthesis/processing chain, which includes the modulators such as LFO's. The only difference is that the bandlimiting must take place (conceptually, this is often not doable in practice) as the final step in an imaginary processing chain (including all modulation etc etc) which happens in CT, the result of which you then sample to form your DT signal.

To take your S&H example: in a CT world, the S&H clock could well sample its input between the samples, and the resulting discontinuous transition would then happen also between the audio samples. If you model your synthesis chain in CT and then bandlimit the whole system, first of all the CT discontinuous transition in the output signal will be slightly smoothed by the anti-alias filter. Further, the sub-sample position of the transition will be reflected in the exact shape of the anti-aliased edge.

Again, whether all of this makes a difference you can actually hear depends on the circumstances, but the point is that a naively sampled LFO does not give you an ideal modulation, but neither does applying modulation from an anti-aliased LFO.


finding period with autocorrelation is not correct I have recorded a signal, which I know is periodic (apart from noise). The period length is unknown. I want to extract the last period from the signal. Before going to a noisy signal, I first tested the method on a self-created noiseless signal. I've read that the autocorrolation is the way to do this, but it did not work for my example. Let me show you what I have done.

First I created a periodic signal of 1000 samples with the Matlab code below:

N_period = 1000;
start_time = 0;
end_time = 1;
dt = (end_time-start_time)/N_period;
fs = 1/dt; %sample frequency
t = start_time:dt:end_time-dt;
y_period = sin(2*pi*t);

I then extended the signal in order to be longer than one period, but smaller than 2 periods:

y = [y_period(0.7*N_period:N_period), y_period];

I know, by construction, that I need to get the last 1000 samples to extract the last period out of the signal.

I first pad the signal with zeros, and then apply the autocorrelation

y_cor = [zeros(1,N), y, zeros(1,N)];
[autocor, lags] = xcorr(y_cor);
figure;
plot(lags, autocor)
grid

As can be seen in the figure below, the second largest peak appears at lag=946 and not at lag=1000.

If the method does not work for a noiseless signal, how can it then work on a noisy signal? Did I do something wrong? How can I improve this?

Actually, implementing my suggestion in the comment leads to a better solution:

N_period = 1000;
start_time = 0;
end_time = 1;
dt = (end_time-start_time)/N_period;
fs = 1/dt; %sample frequency
t = start_time:dt:end_time-dt;
y_period = sin(2*pi*t);

y = [y_period(0.7*N_period:N_period), y_period];

[autocor, lags] = xcorr(y);
figure(1);
plot(lags, autocor ./ (100+length(y) - abs(lags)))
grid

I only calculate the autocorrelation of y without the zero-padding (y_cor) (actually, it does not really matter). As you see, if you divide the autocorrelation by the number of overlapped samples, you get a better result. However, it is not perfect, I've added an (arbitrary) constant 100 to the division to account for boundary effects.


Why can so little digital information be stored on a cassette tape? I had heard that tape is still the best medium for storing large amounts of data. So I figured I can store a relatively large amount of data on a cassette tape. I was thinking of a little project to read/write digital data on a cassette tape from my computer sound card just for the retro feeling. (And perhaps read/write that tape with an Arduino too).
But after reading up about it for a bit it turns out that they can store very small amounts of data. With baud rates varying between 300 to 2400 something between ~200KB to ~1.5MB can be stored on a 90 minute (2x45min) standard cassette tape.
Now I have a lot of problems with understanding why that is.
1- These guys can store 90 minutes of audio. Even if we assume the analog audio quality on them was equivalent of 32Kbps that's about 21MB of data. I have a hard time believing what I listened to was 300bps quality audio.
2- I read about the Kansas City standard and I can't understand why the maximum frequency they're using is 4800Hz yielding a 2400 baud. Tape (according to my internet search) can go up to 15KHz. Why not use 10KHz frequency and achieve higher bauds?
3- Why do all FSK modulations assign a frequency spacing equal to baud rate? In the Kansas example they are using 4800Hz and 2400Hz signals for '1' and '0' bits. In MFSK-16 spacing is equal to baud rate as well.
Why don't they use a MFSK system with a 256-element alphabet? With 20Hz space between each frequency the required bandwidth would be ~5KHZ. We have 10KHz in cassette tape so that should be plenty. Now even if all our symbols were the slowest one (5KHz) we would have 5*8 = 40000 baud. That's 27MB of data. Not too far from the 21MB estimation above.
4- If tape is so bad then how do they store Terabaytes on it?

I had heard that tape is still the best medium for storing large amounts of data.

well, "best" is always a reduction to a single set of optimization parameters (e.g. cost per bit, durability, ...) and isn't ever "universally true".
I can see, for example, that "large" is already a relative term, and for a small office, the optimum solution for backing up "large" amounts of data is a simple hard drive, or a hard drive array.
For a company, backup tapes might be better, depending on how often they need their data back. (Tapes are inherently pretty slow and can't be accessed at "random" points)

So I figured I can store a relatively large amount of data on a cassette tape.

Uh, you might be thinking of a Music Casette, right? Although that's magnetic tape, too, it's definitely not the same tape your first sentence referred to: It's meant to store an analog audio signal with low audible distortion for playback in a least-cost cassette player, not for digital data with low probability of bit error in a computer system.
Also, Music Cassettes are a technology from 1963 (small updates afterwards). Trying to use them for the amounts of data modern computers (even arduinos) deal with sounds like you're complaining your ox cart doesn't do 100 km/h on the autobahn.

But after reading up about it for a bit it turns out that they can store very small amounts of data. With baud rates varying between 300 to 2400 something between ~200KB to ~1.5MB can be stored on a 90 minute (2x45min) standard cassette tape.

Well, so that's a lot of data for when music-cassette-style things were last used with computers (the 1980s).
Also, where do these data rates drop from? That sounds like you're basing your analysis on 1980's technology.

These guys can store 90 minutes of audio. Even if we assume the analog audio quality on them was equivalent of 32Kbps that's about 21MB of data.

32 kb/s of what, exactly? If I play an Opus Voice, Opus Music or MPEG 4 AAC-HE file with a target bitrate of 32 kb/s next to the average audio cassette, I'm not sure the cassette will stand much of a chance, unless you want the "warm audio distortion" that casettes bring – but that's not anything you want to transport digital data.
You must be very careful here, because audio cassette formulations are optimized for specific audio properties. That means your "perceptive" quality has little to do with the "digital data capacity".

I have a hard time believing what I listened to was 300bps quality audio.

again, you're comparing apples to oranges.  Just because someone 40 to 50 years ago wrote a 300 bits per second modem that could reconstruct binary data from audio cassette-stored analog signals, doesn't mean 300 bps is the capacity of the music cassette channel.
That's like saying "my Yorkshire Terrier can run 12 km/h on this racetrack, therefore I can't believe you can't have Formula 1 cars doing 350 km/h on it".

I read about the Kansas City standard and I can't understand why the maximum frequency they're using is 4800Hz yielding a 2400 baud. Tape (according to my internet search) can go up to 15KHz. Why not use 10KHz frequency and achieve higher bauds?

Complexity, and low quality of implementation and tapes. I mean, you're literally trying to argue that what was possible in 1975 is representative for what is possible today. That's 45 years in the past, they didn't come anywhere near theoretical limits.

Why do all FSK modulations assign a frequency spacing equal to baud rate?

They don't. Some do. Most modern FSK modulations don't (they're minimum shift keying standards, instead, where you choose the spacing to be half the symbol rate).

In the Kansas example they are using 4800Hz and 2400Hz signals for '1' and '0' bits. In MFSK-16 spacing is equal to baud rate as well.

Again, 1975 != all things possible today.

Why don't they use a MFSK system with a 256-element alphabet? With 20Hz space between each frequency the required bandwidth would be ~5KHZ. We have 10KHz in cassette tape so that should be plenty. Now even if all our symbols were the slowest one (5KHz) we would have 5*8 = 40000 baud. That's 27MB of data. Not too far from the 21MB estimation above.

Well, it's not that simple, because your system isn't free from  noise and distortion, but as before:

Low cost.

They simply didn't.

If tape is so bad then how do they store Terabaytes on it?

You're comparing completely different types of tapes, and tape drives:
This 100€ LTO-8 data backup tape

vs this cassette tape type, of which child me remembers buying 5-packs at the supermarket for 9.99 DM, which, given retail overhead, probably means the individual cassette was in the < 1 DM range for business customers:

and this 2500€ tape drive stuffed with bleeding edge technology and a metric farkton of error-correction code and other fancy digital technology

vs this 9€ casette thing that is a 1990's least-cost design using components available since the 1970s, which is actually currently being cleared from Conrad's stock because it's so obsolete:

At the end of the 1980s, digital audio became the "obvious next thing", and that was the time the DAT cassette was born, optimized for digital audio storage:

These things, with pretty "old-schooley" technology (by 2020 standards) do 1.3 Gb/s when used as data cassettes (that technology was called DDS but soon parted from the audio recording standards). Anyway, that already totally breaks with the operating principles of the analog audio cassette as you're working with:

in the audio cassette, the read head is fixed, and thus, the bandwidth of the signal is naturally limited by the product of spatial resolution of the magnetic material and the head and the tape speed. There's electronic limits to the first factor, and very mechanical ones to the second (can't shoot a delicate tape at supersonic speeds through a machine standing in your living room that's still affordable, can you).
in DAT, the reading head is mounted on a rotating drum, mounted at a slant to the tape – that way, the speed of the head relative to the tape can be greatly increased, and thus, you get more data onto the same length of tape, at very moderate tape speeds (audio cassete: ~47 mm/s, DAT: ~9 mm/s)
DAT is a digital format by design. This means zero focus was put into making the amplitude response "sound nice despite all imperfections"; instead, extensive error correction was applied (if one is to believe this source, concatenated Reed-Solomon codes of an overall rate of 0.71) and 8b-10b line coding (incorporating further overhead, that should put us at an effective rate of 0.5).

Note how they do line coding on the medium: This is bits-to-tape, directly. Clearly, this leaves room for capacity increases, if one was to use the tape as the analog medium it actually is, and combined that ability with the density-enabling diagonal recording, to use the tape more like an analog noisy channel (and a slightly nonlinear at that) than a perfect 0/1 storage.
Then, you'd not need the 8b-10b line coding. Also, while re-designing the storage, you'd drop the concatenated RS channel code (that's an interesting choice, sadly I couldn't find anything on why they chose to concatenate two RS codes) and directly go for much larger codes – since a tape isn't random access, an LDPC code (a typically 10000s of bits large code) would probably be the modern choice. You'd incorporate neighbor-interference cancellation and pilots to track system changes during playback.
In essence, you'd build something that is closer to a modern hard drive on a different substrate than it would be to an audio cassette; and lo and behold, suddenly you have a very complex device that doesn't resemble your old-timey audio cassette player at all, but a the modern backup tape drive like I've linked to above.


Why the Hilbert transform is not correct using the hilbert function in MATLAB? I am trying to get the phase of a signal with Hilbert transform. However, I found that the unwrapped phase was not correct because the Hilbert transforming using the hilbert function in MATLAB is not correct. The simulation code shown as below:

t=0:0.001:3;
realphase=30*cos(5*t); %signal phase
x=cos(realphase); %signal
H=hilbert(x);

img=imag(H); %Hilbert transform using hilbert function
real=real(H);
ph=unwrap(atan2(img,real)); %unwrapped phase using hilbert function in MATLAB
ht=sin(realphase); %Hilbert transform of cos is sin
pt=unwrap(atan2(ht,real)); %unwrapped phase

subplot(2,1,1)
plot(img)
hold on
plot(ht)
ylabel('Hilbert transform')
legend('imag(H)','sin(realphase)')

subplot(2,1,2)
plot(realphase)
hold on
plot(ph)
plot(pt)
ylabel('phase')
legend('realphase','ph','pt')

The Hilbert transform of $\cos$ should be $\sin$, but with the hilbert function in MATLAB, it is not a $\sin$, why? Could anyone please help me?

I don't think the Hilbert function in Matlab is incorrect - you can establish this by carrying out the same process for a stationary sinusoid (observe the discontinuities at the start/end). Further, the unwrapped phase of the Hilbert transform you've performed does have a sinusoidal characteristic, although with a significant offset.

I think it's more likely you're looking at an interaction of your rapidly changing phase with the windowing implicit to a filter implementation of the Hilbert Transform. Perhaps you could try slowing down the rate in change of phase to see whether you can isolate changes in phase from the window length? You could do this by expanding your time vector (using an increment smaller than 0.001).


Deconvolution of 1D Signals Blurred by a Gaussian Kernel I have convolved a random signal with a a Gaussian and added noise (Poisson noise in this case) to generate a noisy signal. Now I would like to deconvolve this noisy signal to extract the original signal using the same Gaussian.

The problem is that I need a code which does the job of deconvolution in 1D. (I have already found some in 2D but my main aim is 1D).

Can you please suggest me some packages or programs that are able to do so? (Preferably in MATLAB)

Thanks in advance for the help.

I've explained it once on StackOverflow.

Your signal can be represented as a vector, and convolution is multiplication with an N-diagonal matrix (where N is the length of the filter). I am assuming for the sake of the answer that the filter is much smaller than the signal

For example:

Your vector/signal is:

V1
V2
...
Vn

Your filter (convolving element) is:

  [b1 b2 b3];

So the matrix is nxn: (Let it be called A):

[b2 b3 0  0  0  0.... 0]
[b1 b2 b3 0  0  0.... 0]
[0  b1 b2 b3 0  0.... 0]
.....
[0  0  0  0  0  0...b2 b3]

Convolution is:

A*v;

And de-convolution is

A^(-1) * ( A) * v;

Obviously, in some cases de-convolution is not possible. These are the cases when you have singular A. Even matrixes that are not singular, but close to being singular, can be problematic, as they will have large numeric error. You can estimate it by computing the condition number of the matrix.

If A has low condition, you can compute the inverse, and apply it on the result.

Now, let's see some examples in Matlab:

First, I've made a function that computes the convolution matrix.

function A = GetConvolutionMatrix(b,numA)
    A = zeros(numA,numA);
    vec = [b  zeros(1,numA-numel(b))];
    for i=1:size(A,1)
        A(i,:) = circshift(vec,[1 i]);
    end
end

Now, let's try to see what happens with different kernels:

    b = [1 1 1];
    A = GetConvolutionMatrix(b,10);
    disp(cond(A));

The condition number is :

 7.8541

This one is problematic, as expected. After averaging, it is hard to get back the original signal.

Now, let's try some milder averaging:

b = [0.1 0.8 0.1];
A = GetConvolutionMatrix(b,10);
disp(cond(A));

The result is:

1.6667

That goes well with our intuition, mild averaging of original signal is much easier to reverse.

We can also see how the inverse matrix looks like:

 figure;imagesc(inv(A));

Here is one line from the matrix:

  0.0003   -0.0026    0.0208   -0.1640    1.2910   -0.1640    0.0208   -0.0026    0.0003   -0.0001

We can see that most of the energy in each line is concentrated in 3-5 coefficients around the center. Therefore, in order to deconvolve, we can simply convolve the signal again with this approximation:

   [0.0208   -0.1640    1.2910   -0.1640    0.0208]

This kernel looks interesting! It is a sharpening operator. Our intuition is correct, sharpening cancels out blur.


Rotating and Shifting image doesn't change FFT2 I have an image and its fourier transform.
When I rotate it, its fourier transform rotates too, but I can't figure it out. Why does this happen?
On the other hand, when I shift the image, its fourier transform doesn't change. As I know, time shifting means frequency shifting. Am Iwrong?

Time shift corresponds to a phase shift in frequency domain, not a frequency shift. Since you display only the magnitude of the 2D FFT before and after the shift, you do not observe the phase shift.

$$f(x) \Rightarrow F(w) \\ f(x-a) \Rightarrow e^{-j2\pi wa} F(w)
$$

Rotation in space corresponds to a rotation in frequency domain as you observe.

For proofs, see these notes on Some Properties of Fourier Transform


Natural, Forced and Total System Response - Time domain Analysis of DT I am a self-learner of DSP by reading related books.

While reading Time domain Analysis of DT (Discrete Time) system, I suddenly got introduced to difference equation representing a system. It then further stated -

The response of any DT system can be decomposed as :

Total Response = Zero State Response + Zero Input response

Where Zero state response is system response due to input alone and zero input response is when input is zero.

It goes further on to explain Natural, Forced and Total System Response - but unfortunately I am unable to map it to a concept.

For me its like response of a car engine (example per say or anything - a system) is like response when engine is not started + when engine is started - but this does not makes any  sense for me. The engine system why should it consider zero input response at all. Its just a change of its state.

Similarly I am trying to understand Natural / Forced response but do not get the concept quiet clear. Also What is the role of studying them and where can we apply these concepts?

What is meant by recursive and non-recursive system - though my understanding is when system output is fed again into system as its input, it would then make the system call as recursive system - but why does a system need to do that?

First of all, many textbooks use the following definition (but see this answer for a different one):

zero-state response (ZSR) = forced response

zero-input response (ZIR) = natural response

The zero-state response is the response of a system to an input signal given that the initial state of the system is zero. The zero-input response is the output of the system due to a non-zero initial state. For a discrete-time system, the initial state is given by the values stored in memory (i.e. in the delay elements). Since we're dealing with linear systems, the superposition principle can be applied, and the total response of the system is given by ZSR+ZIR. So if the initial state is zero, the total response equals the ZSR, and if there's no input signal, the output is given by the ZIR, which always decays if the system is stable.

In analog circuit theory you have exactly the same. There the initial state is the energy stored in capacitors and inductors. Your car engine analogy is less fortunate because it doesn't really work. You could come up with a mechanical analogy though as long as the system is linear.

It might be a good idea to formulate a new question concerning recursive systems (after having searched the internet and this site). But you're right that a system is recursive if not only the input but also past output values determine the current output. The reason why this is done is the fact that in this way certain system properties can be implemented much more efficiently than by only using the current and past input values. The problem with recursive systems is that they can be unstable, and that quantization effects (in signals and filter coefficients) can more easily result in numerical problems than for non-recursive systems.

PS: If you teach yourself I strongly recommend MIT's open courseware. You should probably start with 6.003 Signals and Systems.


DFT have independent variable in between 0 to $\frac{1}{2}$ no of samples I was learning of Fourier Transform from this book 'Scientist and Engineer Guide to DSP'. Suppose we have a digital periodic signal and there are 128 of samples of it in the time domain i.e X[0] to X[127], then on converting it to frequency domain (only real DFT is considered i.e no complex variables) why do we need only 65 points in the frequency domain, is this related to the Nyquist rate since the highest frequency in the sample can have a maximum value of 64. Are the last 64 points not needed for calculation

This statement is only true if your time-domain signal $x[n]$ is real. A real signal has a Fourier transform that is conjugate-symmetric. Stated differently the "negative frequencies" are just the conjugate of the corresponding "positive frequencies." Therefore, the negative frequencies add no new information, so they can be discarded from the DFT of a real signal.


How do I know if the given signal is discrete? If I am given an expression of a signal then how do I know if it is a discrete signal or not?

If this signal is given how do I know if it's a discrete or continuous?

$$x(n) = e^{j (n \pi/2 + \pi/8) } $$

Fortunately, one typically does not receive signals from strangers in real life. One knows because one is working on an application.

A discrete signal can be equivalent to a continuous signal so you need a set of assumptions tying the two together, like a Nyquist sampling criteria and then you need to test for equivalence under those assumptions.

A function that only exists at discrete points( like time points) and there is no meaningful notion of intermediate values between those points is discrete. Some financial time series are like that.  Arrival times of a Poisson process is another discrete function.  If intermediate values on a continuous range are meaningful, a continuous function is likely.


How to prove energy is preserved in sampling using Parseval's Relation? I have a continuous signal,xc(t) which is bandlimited to =pi/T ,
and x[n] = xc(nT) are samples of this signal.
I need to show there is a Parseval's relation between the original signal and the sampled version.That is I need to prove that energy will be conserved!
Can anyone please help?

Since your signal is band-limited, it has a representation in the orthogonal system of shifted sinc functions,

$$x_c(t)=\sum_{n\in\Bbb Z} x_c(nT)\operatorname{sinc}\left(\tfrac tT-n\right)$$

where $\operatorname{sinc}(x)=\frac{\sin(\pi x)}{\pi x}=\int_{-1/2}^{1/2}e^{i2\pi fx}\,df$.

By the Plancherel identity of the continuous Fourier transform this leads to
$$
\langle \operatorname{sinc}(x-m),\,\operatorname{sinc}(x-n)\rangle
=
\int_{-1/2}^{1/2}e^{i2\pi f(n-m)}\,df
=
\delta_{m,n},
$$
 i.e., orthogonality of the unscaled functions. In the scaled variant above one gets, first by orthogonality and then by rescaling of the time parameter, the Parseval identity for the orthonormal basis $\tfrac1{\sqrt{T}}\operatorname{sinc}(\tfrac tT-n)$, $n\in\Bbb Z$, of the space of band-limited functions,

$$
\|x_c\|_2^2=\sum_{n\in\Bbb Z} |x_c(nT)|^2 \langle \operatorname{sinc}(\tfrac tT-n),\,\operatorname{sinc}(\tfrac tT-n)\rangle=T\,\sum_{n\in\Bbb Z} |x_c(nT)|^2
$$


How does image masking work? I am new to image processing and I am interested how image masking can be used to extract a portion of a image.

Have a look at this OpenCV tutorial.

In the section Application of Mask, they say that you can use a white mask over the area of the image you want to extract.

So I am interested in how this math works. Correct me if I am wrong but if I take the minimum of the pixel from the original pixel and the mask I should get the portion of the portion of the image I want to extract.

This works because the area that I want to extract is covered by the white area of the mask, ie covered by pixels with 255 in all channels, so

min(pixel, 255)

will always give me the pixel value,

while on the other hand the portions I don't want are covered by black pixels(intensity value of 0), so

min(pixel, 0)

will always give me 0 thus I will exclude these sections from my filter.

This my intuition correct ? Is this how image masking works in order to extract the image ?

You may interpret it this way, but this yields a limited interpretation, because for instance the image may have a maximum different from $255$.
I consider masking as a product operation. For a binary mask, the values are $0$ and $1$. Thus if $p$ denotes one pixel value, you get $0\times p = 0$ outside the mask, and $1\times p = p$ inside.

This way, you do not depend on the range of the image anymore. The product can be interpret as the binary AND operation, as you can see in Image masks. The relation between min, max and products could be understood in the more global framework of lattices and Boolean algebras, which are useful in mathematical morphology, and help you define other types of useful masks.

For more continuous operations on "real-valued" images, you can use masks with values in $[0,1]$, whose product with the image allows dimming between white and black mask areas, using intermediate values. For instance in the following, the mask   masks the left part of an image (black mask values), keep the right part, and  affects attenuated values (smooth transition) in the central part.


How to reduce time consuming for calculating pseudo inverse of large matrix in matlab? I have a matrix

P = randn(45875x65536 );
Pi = pinv (P);

I tried to run this code in matlab, but it takes long time

is it possible to split the matrix into smaller matrices then calculate the pseudo inverse?

any suggestion to reduce the time ?

Long story very short: No.

If that was possible, Matlab would likely be doing it already; in any case, the Matlab docs for pinv say that the singular value decomposition method is used; that might not be the fastest method, but it's relatively stable for many cases of P, so it's pretty much desirable unless you know exceptionally well what you're doing or your P has special structure (in which case you'd probably not want to use the Moore-Penrose Pseudoinverse in the first place, maybe?).

Longer story: typically, you multiply your Pi onto something to get the inverse operation to P executed.

To cite Matlab documentation itself:

  You can replace most uses of pinv applied to a vector b, as in pinv(A)*b, with lsqminnorm(A,b) to get the minimum-norm least-squares solution of a system of linear equations. lsqminnorm is generally more efficient than pinv, and it also supports sparse matrices.

But whether you can do that depends on the actual reason you're using the Moore-Penrose Pseudoinverse; not all uses are just replaceable with something else.


Response function for digital camera pixel To help me with doing some signal processing on an image, I'd like to model how a digital camera works.

The camera has a bunch of pixels.  What is the relationship between the value of that pixel in the resulting image, compared to the intensity of the light entering that pixel sensor?

Is it a linear function?
Is there some nonlinear response function?
If it is nonlinear, how  can I recover or model that response function?

If I understand you correctly, what you are asking for is called camera response function (CRF) and in general it is nonlinear and depends on camera device. For fixed device denote its CRF by $f$. Then $f$ maps the set of possible scene irradiances (or illuminances) $\mathcal{I}$ at a given spatial location to the set of possible pixel intensity values $\mathcal{B}$ so $f:\mathcal{I}\rightarrow \mathcal{B}$. In order to estimate $f$ you need to prepare a data set of pairs $(I_j,B_j)_{1\leq j\leq n}$ such that
$$f(I_j) = B_j$$
In practice you can use $I_j$ as raw camera outputs (for instance you may use dng format) of a given scene $S$ and $B_j$ as a pixel values (extracted from jpg image you have taken immediately after dng image) that correspond to them. Then just use any reasonable regression method.

Estimation of CRF is a subject of a large research in computer vision. Let me give you some prominent references. Since I am not expert in this area, I strongly recommend you to consult them.

Debevec, Malik

Grossberg, Nayar


Note segmentation from pitch values My goal is to extract a melodic respresentation for monophonic signal like voice.

I have used YIN algorithm to find the pitch of the signal. Now I want to detect notes from the pitch values (like they will change when a note change and become stable when it does not change). So my goal is to segment pitch values into notes so basically I have to:

split the different sung notes within voiced segments.
assign a constant pitch value to each estimated
note.

So what should be my approach after the pitch estimation ?

(As i have very few DSP experience so I would like a detailed answer. Thank you)

okay, so 9 hours ago i posted in a comment links to other questions/answers regarding this topic for you to maybe pre-digest.

this is the issue of audio-to-MIDI conversion, which has been discussed, particularly among the Music Information Retrieval(MIR) people.  it's not a perfectly solved problem, especially for anything polyphonic.  even picking apart a duet is hard.

for monophonic, it's not so hard to get a decent pitch contour plot (that is the detected pitch vs. time) along with a plot for amplitude vs. time and a plot for pitch confidence or tone harmonicity (or the degree of periodicity of the note waveform) vs. time.  there are other measures that one can do, such as a few different timbre parameters (like brightness, even/odd harmonic energy ratio).

what would be interesting maybe for you to do is, with your favorite monophonic test file (people often use Tom's Diner from Suzanne Vega) and see what pitch contour you get from YIN.  and compare that to the pitch contour you get from using the ASDF (turned into autocorrelation) and peak picking from my other answers:

Fast pitch recognition

How to deal with low fundamental when using AMDF for pitch extraction?

What is an AMDF?

What is the difference between pitch detection and Onset detection?

so first make sure you get a decent pitch contour for your monophonic voice and see if the "octave errors" get patched up before looking for note onsets.

for many notes, the onset is clear by looking at sharp transitions in the amplitude contour curve and in the "harmonicity" curve.  but if you're humming and not having any impulsive or sharp beginnings to notes, these curves make it hard to detect note onsets.  sometimes then you need to look at somewhat sharper pitch transitions and also when the pitch crosses semi-tone boundaries.  but somehow you don't want vibrato to trigger new notes if the vibrato happens to bring your note a little closer to another pitch in the 12 notes per octave scale.

finally, once you do effectively segment your notes, then you need to look at the pitch contour and decide, closer to the note onset, which of the equally-tempered, $\frac1{12}$ octave pitches you are closest to.  but that can have a pitch offset if the whole song is a little sharp or a little flat (but the notes are well tuned to each other).


Public datasets of signals Is there a public database where I can find or anyone can provide me real signals.

I'm looking specifically for pairs of natural signals (e.g. from Sensor networks) to test some cross-correlation and cross-spectrum methods.

The IEEE used to host something called the Signal Processing Information Base, which is gone but there is a mirror at:

http://spib.linse.ufsc.br/database.html

There is also some links from Mike Porter's Ocean Acoustics Library at:

http://oalib.hlsresearch.com/


transfer function of a sampler in the s domain I would like to modelize my whole system into the S-domain. This is a mixed system, there a numerical part (corrector, ADC, DAC) and an analogic part (plant transfer function, sensors, etc...). I know that it is possible to transform a DAC (Digital to Analog converter) or a ZOH into the s domain thanks to the Padé's approximation. Nevertheless I do not know how to do for the ADC as it is non linear like the DAC. Is it possible to do the same with the Padé approximation. Does anyone know the transfer function in the S domain or the approximation of the transfer function of the ADC ?
Thank you and have a nice day !

ADCs don't have a transfer function. An ideal ADC is an ideal sampler with a sampling period T.
Real ADCs are not ideal, they have some sampling delays, jitter, quantization noise, thermal noise, etc. Usually the sampling delays are really small compared to the dynamic of the systems, so we don't model them. Quantization noise and thermal noise should not be a factor unless the signal you measure is really small. And jitter is not an issue unless the signal frequency is really high.
So bottom line, if your ADC is really fast compared to the dynamics of your system and if the measured signal is inside the dynamic range of your ADC, you could pretend that the ADC is not there.


Is there a program to extract features from an audio signal? I am looking for a way to take a .wav file and extract features from it that could be used for further classification (the features could be MFCC statistics, an autoencoder representation, etc.). I know it is likely to be task dependent, but I am looking for some standard set of features that make sense. I am trying to categorize an audio signal (mostly music) into a set of categories such as "quiet music", "thrilling music", etc.

I was wondering if there is existing software (preferably for Linux) that takes such audio signal (such as in a .wav format) and outputs a feature fingerprint.

The openSMILE audio feature extraction toolkit may be able to provide the functionality you desire, where the input is a .wav file and the output extracted audio features. See: http://audeering.com/technology/opensmile/

openSMILE provides a command line executable that is coupled with a configuration file that defines the features to be extracted.  The executable has binaries for Linux and Windows (32 bit), and I've also built the executable from source for use on macOS.  It provides its own objects that calculate features such as MFCCs, energy, pitch, but also "chroma" musical features (which may help with music classification).  It comes with many prebuilt feature configurations that have been used for audio classification tasks, such as the Audio Visual Emotion Challenge (http://sspnet.eu/avec2017/).

Below is an example of how openSMILE can be used to extract features from a single .wav file and a chosen configuration.  The output.arff is a format used by the Weka machine learning library.

SMILExtract -C config/emobase.conf -I input.wav -O output.arff

Custom configurations can also be written using the openSMILE configuration language and also the extensive PDF documentation included with openSMILE's download.

This has the desired functionality as described in the original question:

  looking for a way to take a .wav file and extract features from it that could be used for further classification


Question regarding DTFT of a complex signal I have been doing DTFT practice problems for my DSP course, and I encountered this problem in the textbook that completely stumped me.

The question asks to find the DTFT of the shown signal and to plot its magnitude and phase. I am very comfortable with doing this with real-valued signals, but have no idea where to start when it comes to complex signals. Especially how the plotting would be done. I would appreciate any advice on how to tackle the problem.

I think you almost have the solution. As you suggest in your comment, if you express your signal as :
$$
x(n) = j\cdot \delta(n) - j\cdot \delta(n-4)
$$
As DFT is linear, you will get :
$$
\begin{align}
X(k) &= \sum_{n=0}^{N-1}x(n)\cdot e^{-2j\pi \frac{kn}{N}}\\
 &= j\cdot \sum_{n=0}^{N-1}\delta (n)\cdot e^{-2j\pi \frac{kn}{N}} - j\cdot \sum_{n=0}^{N-1}\delta (n-4)\cdot e^{-2j\pi \frac{kn}{N}}\\
&=j-j\cdot e^{-8j\pi \frac{k}{N}}\\
&=\text{ }...
\end{align}
$$


When does the Amplitude Spectrum equal the Magnitude Spectrum? I had the following question on edX:

I'm failing to understand why the second signal has $M(\omega)=A(\omega)$. First I find the DTFT of the signal:

$$\mathcal{F}\{\delta [n]+\delta [n-1]\}\ =\ 1\ +\ e^{-j\omega}$$

I then deduce the $A(\omega)e^{\ j\phi_A(\omega)}$ representation:

$$1\ +\ e^{-j\omega}\ =\ \big[\cos (\omega)+1\big ]\ -\ j\cdot\sin (\omega)\ =\ 2\cos \bigg (\dfrac{\omega}{2}\bigg )e^{-j\frac{\omega}{2}}$$

The magnitude spectrum is:

$$M(\omega)\ =\ 2\cdot \bigg |\cos \bigg (\dfrac{\omega}{2}\bigg )\bigg |$$

Therefore $M(\omega)\neq A(\omega)$. As a result, $\phi_M(\omega)$ will equal $\phi_A(\omega)$ except contain discontinuities of size $\pm\pi$ at values of $\omega$ where $A(\omega)$ has a zero crossing.

The only scenario I can envisage where $M(\omega)=A(\omega)$ is if $A(\omega)\geq 0$ however as I got the question wrong there must be something wrong with my logic, could somebody clarify?

You are right that also for the second signal you have $M(\omega)\neq A(\omega)$, but the point here probably is that in the fundamental interval $\omega\in (-\pi,\pi)$ the two are identical. It is only at $\omega=\pm\pi$ that the phase $\phi_M(\omega)$ jumps and that $A(\omega)$ changes its sign. So what they probably mean is

$$M(\omega)=A(\omega)\quad\text{and}\quad\phi_M(\omega)=\phi_A(\omega)\quad\bbox[5px,border:2px solid red]{\text{for }\omega\in(-\pi,\pi)}$$

which is correct (and which is not the case for the first signal).

Remember that such quizzes are made by fallible humans.


How to find phase values on precise pixel locations of a given Image? I have computed magnitude and phase spectrum of cameraman Gray scale image  using fft function. Here, we get magnitude and phase spectrum of the whole image. .Please see the reconstruction of image from phase only part from the accepted answer of Image Reconstruction:Phase vs Magnitude  which says  that in an image neighboring pixels of the edge,lines and corner are in phase . So  I want to find phase values on precise pixels

Given an image of dimensions 256*256, how can I find phase values of pixels at locations(100,121),(100,122) and (100,123) ?

Whether it is possible using Fourier Transform or not? If not ,what is the other way to do it? Can anybody explain it with (or without )  code?

clc;
clear all;
close all;

i=imread('C:\Users\RK\Desktop\cameraman.gif');
%i=rgb2gray(i);
i=uint8(i);

figure,
subplot(1,3,1);imshow(i);
title('Cameraman Gray scale Image');

f1=fft2(i);
f2=log(1+fftshift(f1));

m=abs(f2);
subplot(1,3,2);
imshow(m,[]);
title('Magnitude Spectrum');

phase=angle(f2);
subplot(1,3,3);a
imshow(phase,[]);
title('Phase Spectrum');

You can't find the precise phase value of a given pixel because of the uncertainty principle. You can only find the phase value for the pixel's neighborhood using short-time Fourier tranform.

If you want to know the phase for pixel (100,100), ou have to take a neighborhood, for instance (50:150,50:150) and calculate the 2D Fourier transform on this window.


Discrete Fourier Series - Time Index I have been trying to understand Discrete Fourier Time Series (NOT Transform).

It is defined as

$$a_k = \sum_{n=0}^{N-1} x[n]e^{jk\frac{2\pi}{N}n}$$

where N is the Time period and an integer (by definition), and $n$ goes from $0 ... N-1$.

Fair enough. But when N = 0.2 second as the time period of the wave, then how we define $n$, as $N-1$ is negative. I think here we should take $N$ as 200 miliseconds. And then loop $n$ from 0 ... 199.

But then the terms $x[n]$ (the samples) can be defined over either 200 milisec OR $2*10^8$ nanoseconds and so on. So for 200 ms case, we need 200 samples. But lets say I have taken 100 samples distributed equally over 200 ms time period. Now how we define $x[n]$ in this case, should I do like this:  $x[0]=1, x[1]=0, x[2]=1,...$ That is, set alternate values to zero?

Because in the above equation, I need 200 values of $x$.

Or Better (I think):

Define $1 sec = 500ss$, where $ss$ is a hypothetical new unit of time, in that case time period is $N=0.2*500ss=100ss$ (with the condition that the final value is an integer). In this case we can loop over 100 samples easily without worrying over setting alternate values.

Can anyone kindly help me with this or on my thinking above.

Best Regards,

The integer $N$ is the number of samples of $x[n]$ in one period. There need not be any relation to absolute time whatsoever, simply because $x[n]$ are just discrete data that are not necessarily samples of a function of time.

However, if the data $x[n]$ are samples of some continuous-time signal, then the time difference between two adjacent samples is, by definition, the sampling interval $T_s$, which is the inverse of the sampling frequency $f_s=1/T_s$.

In that case, the time corresponding to one period of $x[n]$, i.e., the time difference between two samples $x[n]$ and $x[n+N]$, equals $N\cdot T_s=N/f_s$.


Spectral peak location estimation using complex DFT In this paper a simple method to estimate a spectral peak is proposed, by using quadratic interpolation between three samples of the DFT of the signal. Namely, the position of the peak relative to the index $k$ (they call this difference $\delta$) is calculated with this formula in the paper:
$$\delta = \frac{|X_{k+1}|-|X_{k-1}|}{4|X_{k}|-2|X_{k-1}|-2|X_{k+1}|}$$
This picture will make the formula easier to understand:

Then in the paper it is stated that:

The expression is simple, but it is statistically
biased and performs poorly in the
presence of noise. Some simple changes
[...]
improve its accuracy dramatically, for instance by using the complex
DFT values rather than the magnitudes
as follows:
$$\delta = -\mathrm{Re}\left[\frac{X_{k+1}-X_{k-1}}{2X_{k}-X_{k-1}-X_{k+1}}\right]$$

I understand where the first equation comes from. One can get to that equation by simply finding the coefficients of a parabola such that $y(m) = |X_m|$ for $m\in \{ k-1, k, k+1\}$, and then one can find the maximum for that parabola and get to the first expression for $\delta$. However, I don't understand why when using the complex DFT, the result the paper shows is that one.
Why does that negative sign appear? Why does the denominator seem to be divided by 2? Why do they take the real part instead of the magnitude, completely ignoring the imaginary part?

I admit that I failed my first attempt.
The text suggests that we would take $g(m) = am^2 + b m + c$, interpolating $g(k) = X(k+l)$, for $k \in \{-1, 0, 1\}$, interpolating this would give us polynomial parameterized in $X[k], X[k-1], X[k+1]$, derive, solve and it gives the optimum value at $m = k+\delta$. Not that easy, and it gave me three roots, none of them seemed to be the desired value :(
The comment on the paper is not helpful but they give us a reference [4], from the one of the authors, and there it gives more details. There you will find the equation (3.7) and they give a concise derivation on the Appendix A.
Derivation
Here I will fill some details of the proof.
You express the interpolator polynomial centered at $k + \delta$, so delta will maximize it by construction (more clever the attempt described above).
$$g(m) = a (m - k - \delta)^2 + b$$
Then $g(m)$, determine delta so that interpolates the three samples. Since the quadratic polynomial interpolating three points is unique, we know that $\delta$ here is unique as well.
$$\begin{eqnarray}
X[k-1] &=& g(k-1) &=& b+((k-1) - k - \delta)^2 &=& b + a\cdot(\delta^2+2\delta + 1) \\
X[k  ] &=& g(k  ) &=& b+((k  ) - k - \delta)^2 &=& b + a\cdot(\delta^2) \\
X[k+1] &=& g(k+1) &=& b+((k+1) - k - \delta)^2 &=& b + a\cdot(\delta^2-2\delta + 1)
\end{eqnarray}$$
Start looking at the coefficient of the rightmost terms, and note that
$$2 X[k] - X[k-1] - X[k+1] = 2 a$$
Also note that
$$ X[k-1] - X[k+1] = 4 a \delta$$
$$\delta = \frac{X[k-1] - X[k+1]}{4a} = \frac{X[k-1] - X[k+1]}{2(2X[k] - X[k-1] - X[k+1])}$$
With the sentence

A simple adaptation provides an estimator suitable for estimating tone frequencies from complex DFT outputs.
They introduce the formula we see in the first paper you mentioned.

The adaptation they mention is because $\delta$ is complex, so he takes the real value. Maybe a justification for that would be that $\min_{f} |f - \delta|$ is $Re[\delta]$.


Recording and playback of sound, frequencies question If I record a sound at a rate of 8000 samples per second I know that the highest frequency that can be reproduced from the sound is 4000 Hz. What I want to know is;

  If the sound contains frequencies above 4000 Hz when I play the sound
  back what would I hear and what frequencies would be produced at the
  points where the frequency of the original sound is greater than
  4000Hz.

If your sampling frequency is $f_s=8000$ Hz, your maximum signal frequency is indeed 4000 Hz ($=f_s/2$). If your signal contains frequencies above $=f_s/2$ you would hear the results of aliasing. This means that the original spectrum is folded back into the range $[0,4000]$ Hz. What actually happens is that by sampling with a sampling frequency $f_s$ your spectrum is made periodic with period $f_s$. So if the signal is not band-limited with maximum frequency $f_s/2$, the resulting periodic spectrum will have regions of overlap. E.g., if your signal has frequency components up to 5000 Hz, these frequencies would get folded back to the frequency band between 3000 and 4000 Hz, and consequently lead to distortions in this band. This is the reason why sampling is always preceded by lowpass filtering to attenuate all frequency components that might lead to aliasing (anti-aliasing filter). You can listen to some examples of aliasing here.


Your error might be due to the indices used. Nevertheless, it is always better to state the proper procedure to do it, so that you could find your error.

Let us assume that a 3D scene point $P$ is related to the normalized coordinates of the camera $i$ with the extrinsic parameters $P'_i=R_iP+T_i$ and similarly to camera $j$ with $P'_j=R_jP+T_j$. To find the relative pose we could simply use these two equations to go from the camera coordinates of view $i$ to the real scene and then to the second view $j$. We could then write:

$P=R_i^{-1}(P_i-T_i)$ (from camera $i$ to world)

$P'_j = R_j (R_i^{-1}(P_i-T_i)) +T_j$.  (from world to camera $j$)

If we re-arrange:

$P'_j = R_{rel}P_i +T_{rel}$, where $R_{rel}=R_jR_i^{-1}$ and $T_{rel}=-R_{rel}T_i+T_j$.

where the subscript rel refers to "relative" and is what we seek to find. Keep in mind that, as rotations are orthogonal, $R_i^{-1}=R_i^T$. You can always compute relative orientations in similar manners.

Throughout the computation we used the convention that the transformation $J=[R|T]$ is a mapping from scene to the camera (the inverse is hence the vice versa). The final pose transforms a point in camera $i$ to camera $j$.


Importance of Phase in FFT of an image While processing digital images in the Fourier domain, mostly we exploit the amplitude and not the phase. This could be because the amplitude is much more structured and the amplitude spectrum reveals a lot about the spatial domain as illustrated in HIPR and other answers.
I am really curious to know what does the phase encapsulates, what can we infer from it, what kinds of operations can we perform on the phase of FFT of an image to achieve some objectives?
My specific interest lies in fusing two images in the frequency domain, but any PDF or well-documented links on understanding the phase of FFT(image) will help.
------ Edit (following images as suggested by TimWescott in comments)

The amplitude of an image represents the intensity of the different frequencies in the image. Therefore, it holds the geometrical structure of features in the image (i.e. changes in the spatial domain).
The phase on the other hand, represents the locations of these features (which helps our human eye to better comprehend the image).
Here's a visualization I made in Matlab:

As you can see at the bottom, when we combine amplitude and phase from different images and perform inverse FFT, the resulting image looks much more like the original from which we took the phase. This means that phase holds more information about the image than amplitude.


Correlated signals separation with reference I have a signal S, which needs to be split into two components Sx and Sy.

And I have a signal X, which is a reference signal corresponding to Sx.

I need to perform this split of S and check that resulting Sy is ~Y and Sx ~X (I can use X in the process of filtering\separation, but not Y).

This sounds like a typical denoising task, that can be accomplished with LMS\RLS filters, but in my case signals X and Y are correlated and their BW are overlapping.

They are also non-stationary:
- sometimes X amplitude can decrease to almost 0, then Sx is also 0 and Sy should be estimated as simply S.
- sometimes both X and Y are 0 for a short period of time.
- most of the time X and Y are approx sinusoidal signals with similar central frequencies, one might be slightly shifted w.r.t another.

I tried regular LMS\RLS approaches - assume Sx is noise -> S = Sy + noise, but due to crosscorrelaton between S,X and Y the algorithms best guess is S = Sy.

1) How would you try to solve this?
2) What if we can use Y as well, at least for the beginning? Would it make it simpler?

3) More specific question. Now I have S,X,Y amplitudes in arbitrary units (adc counts). Is it better to scale them? Otherwise, I assume, choice of Sx and Sy will be dependent on amplitude ratios of the signals, or not?

I followed the link from you more current version.  I think your problem is better stated here.

1) I don't think you can solve it as stated.  You will have a one unknown parameter family of solutions.

2) If you can use a portion of Y, and the coefficients remain the same, you can easily solve it.

$$ \vec S = a \vec X + b \vec Y $$

Dot this equation with $ \vec X $ and $ \vec Y $:

$$ ( \vec S \cdot \vec X ) = a ( \vec X \cdot \vec X ) + b ( \vec Y \cdot \vec X ) $$
$$ ( \vec S \cdot \vec Y ) = a ( \vec X \cdot \vec Y ) + b ( \vec Y \cdot \vec Y ) $$

This is a solvable system of two equations, two unknowns.

3) Rescaling will not affect the values of $a$ and $b$.

Hope this helps.

Ced

Followup

If you expect that $a$ and $b$ may be changing, then I think this problem is basically intractable.  At least this approach would be.

These vector equations will apply to any subset of your signal.  You will get the best values for $a$ in sections where you guess that $\vec Y$ is small, and the best values for $b \vec Y$ where $\vec X$ is small (as you already stated).  Also, if you can determine sections where $ \vec X \cdot \vec Y = 0 $ (completely uncorrelated), you will also get good readings on $a$ and $b \vec Y$.  In sections where they are highly correlated, you will get unreliable values for $a$ and $b \vec Y$.

If you have a rough idea of what $a$ is, you can find uncorrelated sections by looking for intervals where $ ( \vec S - a \vec X ) \cdot \vec X \approx 0 $, but this won't get you a better value for $a$ because that's what you started with.

There is no way to separate $b$ from $ b \vec Y $ unless you know something about $ \vec Y $.

Any kind of linear transform is going to look like this:

$$ T( \vec S ) = a T( \vec X ) + b T( \vec Y ) $$

I'm not sure if that could be helpful, I doubt it.


retrieve phase of an extrapolated periodic 1D signal I have a periodic signal on a given grid, say (in matlab format):

t = 1:30;
omega = 2*pi/18.431;
phi = -pi+2*pi*rand(1);   % a random phase [-pi,pi]
x = sin(omega*t+phi);     % the signal
x = x+0.5*rand(1,length(x)); % add some noise

Now I want to retrieve the phase phi. There are a few approaches to that, one is to use fft. The problem is that my grid is not good enough to pick that frequency, hence the phase. (another way is to fit to a sin, but that takes way too long if I need to do it 1e6 times)
I tried to implement the "Autoregressive all-pole model parameters — Burg's method" that is mention in an answer to this question. I manage to pick up the frequency quite well, but fail to pick up the phase.

I was trying to use angle(fft_vec(idx)), where fft_vec is the spectrum obtained after extrapolating the original signal (using arburg), and idx is the index where the frequency omega is picked up. What am I doing wrong? and how can I manage to get that phase by other means?

if you know, in advance, the frequency of the sinusoid, so only the phase (and perhaps the amplitude) is unknown, the simplest noise-immune method i know of to get amplitude and phase is (i'm gonna do this in continuous time and let you translate it to discrete time and eventually to MATLAB):

$$ x(t) = \sin(\omega t + \phi) + n(t) $$

compute

$$ \begin{align}
a & = \int\limits_0^{\frac{2 \pi}{\omega}} x(t) \sin(\omega t) dt \\
& = \int\limits_0^{\frac{2 \pi}{\omega}} \left(\sin(\omega t + \phi) + n(t) \right) \sin(\omega t) dt \\
& = \int\limits_0^{\frac{2 \pi}{\omega}} \sin(\omega t + \phi) \sin(\omega t) dt + \int\limits_0^{\frac{2 \pi}{\omega}} n(t) \sin(\omega t) dt \\
& = \frac{1}{2} \int\limits_0^{\frac{2 \pi}{\omega}} \cos(\phi) dt - \frac{1}{2} \int\limits_0^{\frac{2 \pi}{\omega}} \cos(2\omega t + \phi) dt + \int\limits_0^{\frac{2 \pi}{\omega}} n(t) \sin(\omega t) dt \\
& \approx \frac{1}{2} \int\limits_0^{\frac{2 \pi}{\omega}} \cos(\phi) dt \\
& = \frac{\pi}{\omega} \cos(\phi)  \\
\end{align} $$

and

$$ \begin{align}
b & = \int\limits_0^{\frac{2 \pi}{\omega}} x(t) \cos(\omega t) dt \\
& = \int\limits_0^{\frac{2 \pi}{\omega}} \left(\sin(\omega t + \phi) + n(t) \right) \cos(\omega t) dt \\
& = \int\limits_0^{\frac{2 \pi}{\omega}} \sin(\omega t + \phi) \cos(\omega t) dt + \int\limits_0^{\frac{2 \pi}{\omega}} n(t) \cos(\omega t) dt \\
& = \frac{1}{2} \int\limits_0^{\frac{2 \pi}{\omega}} \sin(\phi) dt + \frac{1}{2} \int\limits_0^{\frac{2 \pi}{\omega}} \sin(2\omega t + \phi) dt + \int\limits_0^{\frac{2 \pi}{\omega}} n(t) \cos(\omega t) dt \\
& \approx \frac{1}{2} \int\limits_0^{\frac{2 \pi}{\omega}} \sin(\phi) dt  \\
& = \frac{\pi}{\omega} \sin(\phi)  \\
\end{align} $$

in both expressions (on the bottom line) we know the middle integral is zero and expect the integral on the right to tend toward zero because $n(t)$ is bipolar and ain't correlated with anything.

i think, then, that you'll find the phase $\phi$ to be:

$$ \phi = \arg\left\{a + jb \right\} $$

which is, if $a>0$

$$ \phi = \arctan\left(\frac{b}{a} \right) $$

otherwise, you might have to add or subtract $\pi$ to that expression.


Multiply a measured quantity of a control-system in the frequency domain during continious operation Currently I am working on a feed-forward control system to actively control vibrations. For that purpose I need to measure accelerations on a vibrating structure and use a transfer function in the frequency domain to calculate a force which should be applied back to the structure.
As the measured acceleration trace is in the time domain I need to perform a fourier transformation, multiply my signal in the frequency domain and do an inverse fourier transform to get back to the time domain to apply the force on the structure using a shaker device.

The whole process looks somewhat like this:

As this whole signal chain should work continously I am unsure how to handle the FFTs? I need a certain signal length in order to perform a fft. So I need to cut my continious signal into chunks and do the ffts of them. Is it usefull to use overlaps with those chunks and how do I determine a suitable FFT length?

Of course I want to get a continious signal at the output. So I need to somehow smooth the cutting points between the chunks after the idft.

So basicaly I know about the individual parts of this process chain, but I have never done such a chain in total and so I am somewhat lost on how to deal with the delay and the stepping I get due to the fft.

Edit:

I think the solution has something to do with the overlap-add method that is talked about in this question:
https://stackoverflow.com/questions/5117839/understanding-overlap-and-add-for-filtering?noredirect=1&lq=1
https://en.wikipedia.org/wiki/Overlap%E2%80%93add_method

I assume that the transfer function that you are using is of a similar form as $V(\eta)$ in the image of your question, with $\eta = \Delta f / f_0$. Often a transfer function is written in terms of the complex Laplace variable $s$, but I assume that your $\eta$ represents some sort of normalized frequency. That is equivalent to using $s = i\,k\,\omega$ with $i$ the unit imaginary number, $k$ some normalization constant and $\omega$ the angular frequency. If you rewrite your transfer function in terms of $s$ instead of $\eta$ you can use one of the many known discretization methods. A simple but good discretization method is the bilinear transform, which requires you to substitute $s$ with the following expression in the transfer function

$$
s = \frac{2}{T} \frac{z - 1}{z + 1}, \tag{1}
$$

with $T$ the sample time of your sensor. After simplifying the resulting (discrete time) transfer function you should have something of the form

$$
V(z) = \frac{b_m\,z^m + b_{m-1}\,z^{m-1} + \cdots + b_0}{z^n + c_{n-1}\,z^{n-1} + \cdots + c_0}, \tag{2}
$$

with $n \geq m$ and all $b_i$ and $c_i$ real scalar constants. When multiplying a discrete signal with the variable $z$ can be interpreted as shifting that signal one sample in time. So given the input signal $a[k]$ (with $k$ an integer used as an index and the time associated with that index could for example be $t = k\,T$) the output $F[k]$ can be related to each other using $(2)$

\begin{align}
(z^n + c_{n-1}\,z^{n-1} + \cdots + c_0)\,F[k] &= (b_m\,z^m + b_{m-1}\,z^{m-1} + \cdots + b_0)\,a[k], \\
F[k+n] + c_{n-1}\,F[k+n-1] + \cdots + c_0\,F[k] &= b_m\,a[k+m] + b_{m-1}\,a[k+m-1] + \cdots + b_0\,a[k].
\end{align}

Solving this for $F[k]$ gives

$$
F[k] = b_m\,a[k+m-n] + b_{m-1}\,a[k+m-1-n] + \cdots + b_0\,a[k-n] - c_{n-1}\,F[k-1] - \cdots - c_0\,F[k-n],
$$

so a weighted sum of previous inputs and outputs. This relation for $F[k]$ allows you to calculate a new output every time you get a new input, so no need to calculate it using chunks. I can be noted that initially you do not know the past values of the input and output. However, usually one set those to zero or some other initial value, which should only change the initial transient behavior of the output (assuming that the filter is stable).


Why Z-transform is considered as separate transform? The mathematical formula of the Laplace and Z transforms are same with just one difference.  I.e. in the first we use $t$ for continuous-time signal and in the latter uses $n$ for discrete-time signal. I don't think that there are any other differences.

While discussing Fourier transform for continuous-time (normally analog) signals, we use the continuous Fourier transform and for discrete-time (most often digital) signals we use the discrete-time Fourier transform (DTFT).

So my question is why the Z transform is considered a separate transform? Why it is not named the "Discrete Laplace Transform" ?

There is indeed a transform called discrete Laplace transform and it is of course closely related to the $\mathcal{Z}$-transform. The (unilateral) discrete Laplace transform of a sequence $f_n$ is defined by (cf. link)

$$\mathcal{L}_T\{f_n\}=\sum_{n=0}^{\infty}f_ne^{-snT}$$

with some $T>0$. The discrete Laplace transform can be interpreted as the Laplace transform of a sampled function $f(t)\cdot\sum_n\delta(t-nT)$  with $f_n=f(nT)$.

In practice it is not convenient to have the factor $e^{sT}$ appear everywhere. If one substitutes $z=e^{sT}$, the discrete Laplace transform is called (unilateral) $\mathcal{Z}$-transform:

$$\mathcal{Z}\{f_n\}=\sum_{n=0}^{\infty}f_nz^{-n}$$

The same can obviously be done for the bilateral versions of the transforms, where the integrals and the sums start at $-\infty$.


Does a phase rotation introduce a delay? Suppose I have 2 complex time series, $x_1$ and $x_2$.

If I measure the delay between $x_1$ and $x_2$, I get $t$.

If I make $x_2 = x_2 e^{-j 2 \pi 10}$, and measure the delay between $x_1$ and $x_2$, do I still get $t$?

I would assume a discrete-time context (though your notation suggests continuous-time).

Then you have a misunderstanding of the relation between the "delay" and "phase". Delay is measured in time-domain whereas the phase rotation is defined in frequency-domain.

So if there was a delay $\pm d$ between $x_1[n]$ and $x_2[n]$, then $y_1[n]$ (defined as the inverse DTFT of $Y_1(\omega)$ defined below) will have zero delay wrt $x_2[n]$ in time:

$$ Y_1(\omega) = e^{ \mp j \omega d } X_1(\omega)$$


image processing vs Computer vision? I have already tried searching about subject and the only one thing I could understand is that computer vision involves object recognition but image processing does not involve object recognition

What are other differences between image processing & computer vision??

Digital image processing is an extension of digital signal processing and linear system theory into two dimensional signals.

Image processing involves all low level tasks such as filter design and filtering, spatial scaling, sampling, intensity manipulations, geometry manipulations, Fourier analysis and spectrum analysis, motion estimation, noise reduction, edge detection, image enhancement and restorations etc.

Computer vision involves high levels tasks such as object recognition, scene recognition and image understanding; a more computer science perspective than a DSP type. Note that computer vision inevitably will use image processing at its preprocessing or post processing steps.


Relationship between Carrier Frequency Offset and the Constellation Diagram As part of a project I am doing, I have to investigate how having a carrier frequency offset would affect the constellation pattern of the received signal.This question has got me baffled because to me a carrier frequency offset should not affect the constellation pattern because any offset would be dealt with on the receiver side.I might not be understanding some concepts very clearly.

ie.  On Transmit side we transmit

$$
    s(t) = x(t)\cos(2 \pi (f_c + f _\phi) t)
$$

then on receive side

$$
    r(t) = 2 \text { LPF  } (s(t) \cos(2 \pi (f_c + f_\phi)) t) = x(t)
$$

So regardless of the offset we have on the carrier frequency there should no be no effect on the constellation diagram.Any help will be appreciated?

I think that you misunderstand the problem. The actual problem is that the frequency offset is unknown to the receiver. I won't give away the answer since this appears to be a homework project, but I'll try to help you by giving a hint. First, imagine what would happen if there was an unknown carrier phase offset at the receiver. How would this affect the constellation? If you know the answer to that question, then think of a carrier frequency offset as a continuing change in phase offset. So how would this affect the constellation?


PLL: cycle slipping correction i am working on a receiver in the space communication. Under low SNR, big frequency offset, phase offset. I started with phase/frequency correction. I used digital PLL and it doesn't work for it. I have read in the book by Rice about "cycle slip" and I can use slip detector in such case. But I dont understand how such detector can be implemented. Can someone provide some information /explanation about it?
Does anyone know the simplest algorithm for cycle slipping correction?
I have found in the book by Rice ( Appendix C) about the cycle slip -effect. I have never heard about it, but it describes perfectly what I get in my simulation of phase correction. Rice didn't provide information about how to correct this effect but mentioned it.

Costas loop (and other PLL-based synchronization algorithms) tend to lose lock when the frequency offset is higher than the pull-in range (i.e. maximum frequency offset the said algorithm can cope with). Thus, a coarse frequency recovery algorithm is typically included in practice. The general practice is to put the coarse frequency recovery before the Costas loop. There are many ways of doing coarse frequency recovery including FFT-based algorithms and the frequency-locked loop (FLL) (described in Harris paper. My favorite is FFT due to its generality. The FFT algorithm can be found on the GOES satellite user manual and is summarized in the figure below.

The algorithm works by:

Carrier detection by absolute-law (residual carrier), square-law (BPSK) or fourth-power law (QPSK). Square law detector involves squaring the IQ signal while fourth-power law uses the fourth-power of the IQ signal.
Finding the FFT
Averaging successive FFT sequences
Finding the peak by looking at 8 adjacent bins with the highest magnitude.

Below is a GNU Radio example for BPSK (square law). As you can see, an FFT is included at the very beginning to reduce the frequency offset (in my case it was Doppler shift/rate + carrier oscillator offset).

If you are interested, you could have a look at a paper I wrote a couple of years ago. The paper covers all aspects of synchronization including phase, frequency, and timing.
I hope this helps. Good luck.


Bridging CTFT and DTFT for a cosine I'm trying to understand how I can start from the CTFT of a signal and end up with a DTFT.

For example if I take a basic example:

$$\begin{aligned}
x(t) &= \cos(\omega_x \cdot t) = \frac{1}{2} \cdot \left( e^{j\omega_x t} + e^{-j\omega_x t} \right) \\ \implies X(\omega) &= \pi \cdot (\delta(\omega - \omega_x) + \delta(\omega+\omega_x)) \\ x_c(t) &= x(t) \cdot \sum_{n=-\infty}^{\infty}\delta(t - nT_s) = \sum_{n=-\infty}^{\infty}x(nT_s)\delta(t-nT_s) \\ \implies X_c(\omega) &= X(\omega) * \left( \omega_s \sum_{n=-\infty}^{\infty} \delta(\omega - n\omega_s)\right) \\
\end{aligned}$$

$$ X_c(\omega) = \omega_s \pi \sum_{n=-\infty}^{\infty} \left(  \delta(\omega - \omega_x - n\omega_s) +  \delta(\omega + \omega_x - n\omega_s) \right)
\tag{1} \label{1} $$

From there I'm lost and everything crumbles. I'm only trying to get the DTFT of a cosine which is:
$$
\cos(\Omega_0 n) \Longleftrightarrow \pi \sum_{n=-\infty}^{\infty} \left(  \delta(\Omega - \Omega_0 - n2\pi) +  \delta(\Omega + \Omega_0 - n2\pi) \right)
\tag{2}
\label{2}
$$

How can I obtain $\eqref{2}$ starting from $\eqref{1}$?

I hope what I'm trying even makes sense. After all the DTFT with infinite period is the CTFT so I suppose there's a link we can make between these equations?

Thanks

You're almost there. First, there's a small scaling error in the transform $X_c(\omega)$. Since the Fourier transform of an impulse train is given by

$$\mathcal{F}\left\{\sum_{n=-\infty}^{\infty}\delta(t-nT)\right\}=\omega_s\sum_{k=-\infty}^{\infty}\delta(\omega-k\omega_s)\tag{1}$$

with $\omega_s=2\pi/T$, and since multiplication in the time domain corresponds to convolution in the frequency domain with a scaling of $1/2\pi$ (if you use $\omega$ as the independent frequency variable), the transform of the sampled signal is

$$X_c(\omega) = \frac{\omega_s}{2} \sum_{n=-\infty}^{\infty} \big[  \delta(\omega - \omega_x - n\omega_s) +  \delta(\omega + \omega_x - n\omega_s) \big]\tag{2}$$

Now you just need to use the relation between $\omega$ and $\Omega$

$$\omega=\frac{\omega_s}{2\pi}\Omega\tag{3}$$

and the scaling property of the Dirac Delta impulse

$$\delta(ax)=\frac{1}{|a|}\delta(x)\tag{4}$$

which for your example results in

$$\delta(\omega)=\frac{2\pi}{\omega_s}\delta(\Omega)\tag{5}$$

and you'll arrive at the correct expression for the DTFT of a sampled cosine.

Also take a look at this answer to a related question.


Can't implement correctly Python IIR and FIR Filters I'm trying to implement 2 lowpass filters and see the difference between them from plots( one FIR filter that needs to be designed with the Hann window method and another IIR filter using Chebyshev type I). The problem is that I can't make 2 correct plots to see the difference between them and I don't know exactly if they are correctly implemented...I have fc1,fc2, ft, rp and rs as parameters. Here it's my python code:

 import matplotlib.pyplot as plt
 import scipy.signal as signal
 import numpy as np

def plotInTime(x,fs):
    t=np.arange(0,np.size(x,0)/fs,1/fs)
    plt.figure()
    plt.plot(t,x)
    plt.xlabel('Time (s)')
    plt.ylabel('Amplitude')
    plt.grid(True)

def plotInFrequency(x,fs):
    N=int(np.size(x,0)/2)
    if np.size(x,0)==1:
        X=np.fft.fft(x,axis=1)
    else:
        X=np.fft.fft(x,axis=0)
    X=np.abs(X)
    X=X[:N]
    f=np.arange(0,fs/2,fs/2/N)
    plt.figure()
    plt.plot(f,X)
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Magnitude')
    plt.grid(True)

fs, x = readWav('Ring05c.wav')

plt.close('all')
plt.figure(figsize=[12,12])
plt.subplot(3,2,1)
plotInTime(x,fs)
plt.title('Original Signal')
plt.subplot(3,2,2)
myDSP.plotInFrequency(x,fs)
plt.title('Spectrum for the Original Signal')

fc1=1200
fc2=3000
fc= np.arange(fc1,fc2+1)
ft=150
rp=0.1
rs=40
t = np.linspace(-1, 1, 201)
L, Wn = signal.cheb1ord(fc1/(fs/2), fc2/(fs/2), rp, rs);
b, a = signal.cheby1(L, 1, Wn,'low')
h = signal.lfilter(b, a,x,axis=0)
plotInFrequnecy(h,fs)
w=signal.hann(65) #idk what to do after this

Right now I have an error that says: Error: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

    import myDSP
        import matplotlib.pyplot as plt
        import scipy.signal as signal
        import numpy as np

        fs, x = myDSP.readWav('Ring05c.wav')

        fc1=1200
        fc2=3000
        fc= np.arange(fc1,fc2+1)
        ft=150
        rp=0.1
        rs=40
        t = np.linspace(-1, 1, 201)
        #x = (np.sin(2*np.pi*0.75*t*(1-t) + 2.1)+0.1*np.sin(2*np.pi*1.25*t + 1) +0.18*np.cos(2*np.pi*3.85*t))
wp=np.array([fc-ft/2])/(fs/2)
ws=np.array([fc+ft/2])/(fs/2)
        L, Wn = signal.cheb1ord(wp,ws, rp, rs);
        b, a = signal.cheby1(L, rp, Wn,'low')
        f,H=signal.freqz(b,a,worN=256,plot=None,fs=fs)
        plt.plot(b,'.-b',label='b coefficients')
        plt.plot(a,'.-r',label='a coefficients')
        plt.legend()

        #h = signal.lfilter(b, a,x,axis=0)
        #w=signal.hann(65)