real time signal analysis of output from an audio device I'm writing a test program that can detect pops/cracks/noise (unknown frequencies and duration) in the audio output of an audio device. The purpose of this program is to reduce human error and make the test more systematic.

The test consists of playing a pure sine tone (say, 1kHz) by the test device, and recording the analog output with a recording device. The minimum test duration is 8 hours and any artifacts present in the output will be logged and saved for future access.

Currently, my approach is as follows:

Sample Rate(playback device): 48000Hz

Sample Rate(recording device): 48000Hz

Windowing function : Hanning window

FFT size = total samples

Record 10 seconds of audio samples.
Apply window, run FFT on the recorded samples.
Record the next 10 seconds of samples.
Apply window, run FFT on new samples.
Compare both FFTs. If the frequency bins have similar power (to a set episilon), then the audio signal is clean. Otherwise, log the timestamp of this occurence.
Repeat steps 1-5 until 8 hours have surpassed.

A 1kHz sine tone, along with a 0.01s 3.5kHz sine tone that repeats randomly, is used as a test signal for the program. Unfortunately, my results have a high frequency resolution and it seems difficult to detect what is considered noise.

I am not well-versed in DSP and was wondering how I can improve the accuracy of my results?

I'm writing a test program that can detect pops/cracks/noise (unknown frequencies and duration) in the audio output of an audio device.

Pops/cracks are impulsive in nature and therefore quite broadband. noise is very broad but I would assume here that is anything additive to the sinusoid.

You can catch all three, given reasonable assumptions, using the same way that Total Harmonic Distortion (THD) is measured.

To do that, setup a very narrow band-stop filter, also known as a Notch filter, and apply it to the recorded signal directly. The cut-off frequency of the Notch filter should be exactly on the frequency of the sinusoid you are recording.

The signal at the output of the notch filter will contain everything except the sinusoid signal.

Which means that you would be able to capture cracks and pops, because they would spill across the spectrum as well as any additive noise to the signal.

To quantify the extent to which your signal is "noisy", you can then apply an envelope detector, which in your case is basically a moving average filter followed by a threshold.

The idea there is that you integrate the "residual signal" (what is left after you have removed the sinusoid) over a period of time depending on the sort of temporal resolution you want to achieve and then you threshold that quantity. Basically, if the total sum of the signal over 5 seconds (for instance) becomes larger than some threshold (similar to the $\epsilon$ you seem to already have used), then that time interval is marked as "noisy".

For a notch filter you can design one using plain simple pole-zero placement or use a Twin-T.

If the Q-factor of those filters is not satisfactory, in other words, if the Notch filters you obtain with these methods are not selective enough, then you can simply subtract the sinusoid at the input from the recorded signal. Of course, this has to be done "carefully" because any phase difference inserted in the middle would result in an offset of your $\epsilon$ and to cut a long story short, prior to subtracting the sinusoid, you would first have to make sure that its phase is aligned, in which case we are into the Phase Locked Loop land. The PLL will basically be generating a local sinusoid that is locked on to the phase and frequency of the incoming sinusoid and you can use that to subtract it from your recorded signal. This would make it probably the sharpest of notch filters.

You could of course sum all the bins except the bin that contains the sinusoid you are trying to measure but the challenging aspect there is going to be the resolution of the DFT.

Hope this helps.

EDIT:

Just as a note, I thought I would just say here that what you are doing there should be done in a multi-threaded application, so that the recording thread proceeds uninterrupted. When you write:

  Record 10 seconds of audio samples.
  Apply window, run FFT on the recorded samples.
  Record the next 10 seconds of samples.

then, it is good to keep in mind that a single thread application would miss part of "reality" between steps 1 and 3 when it was processing the frame (but not recording). You can minimise this "gap" error by processing smaller chunks of audio but within limits. To truly do it in as close to real-time as possible way, you would need to have one thread recording sound and stuffing a queue with frames and another thread on the other side of the queue popping frames for processing.


Maximum Likelihood Detection of Signal Vectors in Gaussian Noise Consider a binary-input additive white Gaussian noise channel.

Let $\mathbf{x}_0 = (\sqrt{Es},\sqrt{E_s},⋯,\sqrt{E_s})$ and $\mathbf{x}_1 = (-\sqrt{E_s},-\sqrt{E_s},⋯,-\sqrt{E_s})$, be two codewords of length $d$. Suppose $\mathbf{x}_0$ is transmitted and the received vector is given by $\mathbf{y} = \mathbf{x}_0 + \mathbf{n}$, where the noise vector $\mathbf{n}$ is an i.i.d. Gaussian random vector with  mean zero and variance $\frac{N_0}{2}$ . Suppose that maximum-likelihood decoding is employed. The pairwise error probability $P_d$ is the probability that the decoder chooses the incorrect codeword $\mathbf{x}_1$ as the decision instead of the originally transmitted $\mathbf{x}_0$. Show that  $P_d=Q(\sqrt{\frac{2dE_s}{N_0}})$, where $Q(x) = (\frac{1}{\sqrt{2 \pi}})\int ^{\infty}_{x}e^{-\frac{t^2}{2}}dt$

I know how to calculate the error probability for bpsk in the other type,but not codewords type,does anyone know how to calculate the probability for codeword type?

The probability of error of the ML detector is equal to the probability that the received vector $\mathbf{y}$ is closer to $\mathbf{x}_1$ than to $\mathbf{x}_0$, which is equal to the probability that the noise component in the direction of $\mathbf{x}_0-\mathbf{x}_1$ is greater than half of the Euclidean distance between $\mathbf{x}_0$ and $\mathbf{x}_1$. The Euclidean distance between $\mathbf{x}_0$ and $\mathbf{x}_1$ is

$$D=||\mathbf{x}_0-\mathbf{x}_1||=\sqrt{\sum_{i=1}^d(x_{0,i}-x_{1,i})^2}=\sqrt{\sum_{i=1}^d4E_s}=2\sqrt{dE_s}\tag{1}$$

The noise variance in the direction of $\mathbf{x}_0-\mathbf{x}_1$ equals $N_0/2$ (as in any other direction). The probability that a zero mean Gaussian noise variable with variance $N_0/2$ assumes a value greater than $D/2$ is given by

$$P_E=Q\left(\frac{D/2}{\sqrt{N_0/2}}\right)=Q\left(\frac{\sqrt{dE_s}}{\sqrt{N_0/2}}\right)=Q\left(\sqrt{\frac{2dE_s}{N_0}}\right)\tag{2}$$


Relationship between "received" signal power and "transmitted" energy per bit? In digital communications textbooks (e.g., Digital Communications by Bernard Sklar), the relationship between the "received" signal power $P_r$ and the "transmitted" energy per bit $E_b$ is given as follows: $P_r=E_b R_b$ where $R_b$ is the bit rate. What confuses me in this equation is that the channel attenuation is not considered in the equation. Shouldn't the equation be as  $P_r=\alpha^2E_b R_b$ instead where $\alpha$ is the channel attenuation?

You are correct that it could be argued that every physical channel will introduce attenuation. In digital communications, however, it is common to define "channel models" that simplify the equations and make the problem easier to understand and to solve.

The most popular model is the "additve white Gaussian noise" (AWGN) model, where the only effect of the channel is to introduce noise. In other words, there is no attenuation. In that model, Sklar's relationship between power and bit energy is correct. This model is very popular because, in many cases, you can ignore the attenuation by assuming an appropriate gain mechanism in the receiver.

Of course, in many other cases the AWGN model is not appropriate. In the wireless channel, for instance, the channel gain is so fundamental that it gets incorporated into the model. If the channel is narrow-band, then the relationship between power and bit energy is similar to what you propose, with $\alpha$ random and possible time-varying.


Detect presence of a specific frequency (or lack thereof) I know very little DSP beyond an intro undergrad course. I was looking to write a C algorithm that takes ADC samples and determines if a specific frequency is present (in the range of 100KHz); not caring about any other components of the signal.

I'm not sure where to start. I looked at Goertzel algorithm; would this be a recommendable approach?

I was wondering if it is feasible to analyse a signal at 100KHz with a ARM micro controller, or will I need a special DSP breed?

Should I even be tackling this with a PLL or something similar? If so, does anyone have some design references?

Thanks!

Regarding the implementation question, maybe you could tell us more about your hardware: ARM family (cortex M3 or M4?), clock rate...

A fixed-point implementation of the Goertzel algorithm is a dozen instructions, so this would indeed be doable with a sample rate of 1 MHz on a low-end ARM micro (say STM32F1 or LPC1343) clocked at 72 MHz (assuming you don't have much overhead for reading your ADC, eg. by letting the DMA subsystem handle this).

However:

You might hit the upper range of your ADC sample rate, and/or will have little margin to design your anti-aliasing filter (say 105kHz is the upper frequency of interest... You can sample at 210kHz but will need the steepest anti-aliasing filter... Or at 800kHz and have more margin for the transition band of your anti-aliasing filter).
This still looks like a waste of CPU cycles, especially if you are only interested in a narrow band of the signal. Keep in mind that if the only thing that interests you is the 99kHz to 101kHz band, you actually need a sample rate of just 4kHz to capture all the useful information - provided you have a mean of modulating the input signal by a 100kHz sine wave.

In such situations, and if you control the target hardware, a bit of pre-processing in the analog domain could be handy, such as:

Frequency shifting your signal by 100kHz so the band of interest becomes the base-band and so that you can sample it at a low sample-rate (requires a sine oscillator, a multiplier, and a low-pass filter).
Applying a narrow band-pass filter in the area of interest + rectifier + low-pass filter - this would measure the energy of the frequency band of interest.


Indexing in DFT (from an old paper) There is a nice paper on explaining DFT from the 1960s in IEEE A guided tour of the fast Fourier transform. The author uses the following definitions of DFT
DFT
$$
X(j)=\sum_{k=0}^{N-1} x(k) \exp \left(-i 2 \pi\left(\frac{j}{N}\right) k\right)
$$
Inverse
$$
x(k)=\frac{1}{N} \sum_{j=0}^{N-1} X(j) \exp \left(i 2 \pi\left(\frac{j}{N}\right) k\right)
$$
where the indexes j = 0, 1, 2, ..., N-1 and similarly k=0, 1, 2,..., N-1.
Now the authors show a figure, where the j and k indices run from 0 to N not N-1. Let us say we had 10 data points, so N=10; and j and k should run from 0 to 9 not 10. Is this a typographical error in the figure?
It seems that his N also starts from zero, then the figure is consistent but the summation formula has N-1.

The figures are correct. You can see that there are samples at indices $0,1,\ldots,N-1$. There are no samples shown at index $N$, neither in the time domain nor in the frequency domain. The value $N$ is only shown on the abscissa because in the frequency domain it corresponds to the sampling frequency, and in the time domain it corresponds to the end of the continuous-time signal represented by $N$ samples (if we assume that each sample represents a portion of length $\Delta T$).
For a very detailed discussion on the definition of the actual duration of a discrete-time sequence have a look at the answers to this question.


Kernel Convolution in Frequency Domain - Cyclic Padding I don't know whether this is the right place to post this, but I suppose it is.
I know that frequency multiplication = circular convolution in time space for discrete signals (vectors).
I also know that "the convolution theorem yields the desired linear convolution result only if $x(n)$ and $h(n)$ are padded with zeros prior to the DFT such that their respective lengths are $N_x+N_h-1$, essentially zeroing out all circular artifacts."
and everything works with vectors.. but my goal is circular convolution with matrices as in this paper:
Victor Podlozhnyuk (nVidia) - FFT Based 2D Convolution
If you watch the first two figures (figure 1 and 2) you'll see that the kernel is padded in a weird way I've never seen before, what's this?

In my StackExchange Signal Processing Q38542 GitHub Repository (See Applying Image Filtering (Circular Convolution) in Frequency Domain) you will be able to see a code which implements 2D Circular Convolution both in Spatial and Frequency Domain.

Pay attention to the function CircularExtension2D().
This function align the axis origin between the image and the kernel before working in the Frequency Domain.
Remember that for Discrete Signals the implicit assumption on signals, In frequency Domain analysis, is being periodic (Circular).
This is basically what's in FFT Based 2D Cyclic Convolution.

In the discrete case one could indeed apply Circular Convolution by element wise multiplication in the Frequency Domain.

With proper padding one could apply linear convolution using circular convolution hence Linear Convolution can also be achieved using multiplication in the Frequency Domain.

See:

Regarding your questions:

The filter is just an array of numbers. As long as you are after 2D Circular Convolution there is no constraints on the Filter. If it is valid for 2D Spatial Circular Convolution it is valid for Frequency Domain Circular Convolution.
Up to Floating Point Quantization errors both are mathematically equivalent (See Convolution Theorem).
If the Convolution Kernel is similar in its size to the image and both are large enough Frequency Domain Convolution becomes more efficient than Spatial Domain.


Isolate the non blurred part of foucsed image I have an image taken by a mobile phone camera with focus on a part of it let say a face or whatever.
I want to know -approximately-  where is the focused area. For example, a bounding box around the focused area.
example :

In general, this is an image segmentation problem (http://en.wikipedia.org/wiki/Image_segmentation) into which you would be trying to isolate the focused to the non-focused regions of the image.

Optical lenses are equivalent to low pass filters anyway and the effect of a low pass filter on a signal is to smooth it out by limiting the higher frequency content.

In digital image processing, low frequencies correspond to big objects (general form) and higher frequencies correspond to small objects (details).

Therefore, to isolate the focused region, you could be looking for areas that are at a higher standard deviation relative to the blurred areas which would be at a lower standard deviation because of the smoothing operation of the out-of-focus lens.

To distinguish 'blurred' vs 'non-blurred' regions  you can setup a "window" of dimensions $M \times N$ and slide it over your image.

At each position of the window compute the variance of the "patch" values as if they were a one-dimensional vector of $M \times N$ values. This will result in another "image" (or more generally a two-dimensional array) depicting the variation of the standard deviation at each location on your image. You can then collect all of these values and obtain their histogram. The histogram will exhibit a multimodal distribution (For more information please see http://en.wikipedia.org/wiki/Multimodal_distribution ).

In the case of the image in your post, the histogram is likely to exhibit two modes. One corresponding to relatively high values of standard deviation (focused area) and one corresponding to relatively low values of standard deviation (blurred areas).

You can then position a threshold, for example just between the two modes to establish a criterion between your focused and unfocused regions and then select all "patches" (or MxN window areas) whose standard deviation is higher than your threshold. (For more information please see: http://en.wikipedia.org/wiki/Thresholding_%28image_processing%29)

This will isolate the focused parts of the image for you.

Hope this helps.


Reconstructing a sine wave from an interval shorter than half its wavelength I know that the function that I search is a sine wave of the form
$$A \cdot \sin(x) $$
where my $A$ is unknown. I have samples of the sine wave at discrete points within an interval that's shorter than $\lambda/2$ represented by the red line in the plot.

It must be possible to somehow determine $A=1$ through reconstruction or so from my samples given by the red line. How can i do that?

If your signal is really as simple as
$$x(t)=A\sin(\omega_0t)\tag{1}$$
with known $\omega_0$, and you have observations $y(t_i)$, which are noisy samples of $x(t)$ at known time instances $t_i$, then a simple solution would be the least squares estimate
$$\hat{A}=\frac{\displaystyle\sum_iy(t_i)\sin(\omega_0t_i)}{\displaystyle\sum_i\sin^2(\omega_0t_i)}\tag{2}$$
Of course, this simple solution won't work if your signal actually has the form
$$x(t)=A\sin(\omega_0t+\phi)+c \tag{3}$$
with unknown phase $\phi$ and DC-offset $c$. However, you can compute optimum least squares estimates for that problem too. This is discussed in Cedron's answer.


What are the parameters of transducer that I can get from an Anechoic chamber? I have a bunch of piezo elements, that I want to use as "Mics", for that I will test their reaction in a Anechoic chamber. My question is what the are characteristics that I can measure that describe those elements and can help me analysing the generated sound later, I already have in mind the following :

Amplitude response
Phase response

Thanks for any help in advance !

Mainly frequency response as a function of the 3D direction of incidence. Also known as "polar pattern" or "polar balloon". You can also look at distortion and noise but that's more tricky.


Carrier Frequency Offset, Rotating Constellation and Correlator I have been further reading up on carrier recovery methods and something is not clear to me:

It's written that a carrier frequency offset results in a rotating constellation diagram but if you are sampling at the output of a correlator such as RRC it's not clear to my why this is the case.

For an all digital receiver, the ADC samples will be rotated prior to going into the correlator but it appears the correlation process in conjunction with the frequency offset will do more than rotate the constellation diagram.

Am I wrong?

The correlator is just a filtering operation. So in the case of passband pulse-amplitude modulation (PAM), such as QPSK or QAM, the (noiseless) received analytic signal before demodulation and sampling can be written as
$$r(t)=e^{j((\omega_c+\Delta\omega) t+\phi))}\sum_{k}A_kp(t-kT)\tag{1}$$
where $\Delta\omega$ is a frequency offset, $\phi$ is a phase offset, $\omega_c$ is the carrier frequency as estimated by the receiver, $A_k$ are the complex data symbols, and $T$ is the symbol period. The pulse $p(t)$ represents the combination of the transmit filter, linear distortions on the channel, and the receive filter (correlator).
Demodulation with $e^{-j\omega_ct}$ and sampling at the symbol rate $1/T$ gives
$$\tilde{r}[n]=e^{j(\Delta\omega nT +\phi)}\sum_kA_kp((n-k)T)\tag{2}$$
If for simplicity we assume that the combined pulse $p(t)$ satisfies the Nyquist criterion $p(kT)=\delta[k]$ (i.e., there is no inter-symbol interference), then $(2)$ simplifies to
$$\tilde{r}[n]=e^{j(\Delta\omega nT +\phi)}A_n\tag{3}$$
which should make it clear that at time $t=nT$ the received symbol is a rotated version of the original symbol, where the rotation has a constant component $\phi$ (the phase offset), and a time-dependent component $\Delta\omega nT$ caused by the frequency offset. Consequently, if $\Delta\omega\neq 0$, the constellation rotates.


Doubling of the scipy.signal's bandpass filter order during the design I am trying to extract filter coefficients for a bandpass, using the "scipy.signal.iirfilter" function:
b, a = filt = iirfilter(
        N=2,
        btype='bandpass',
        ftype='cheby2',
        Wn=np.array([3000, 4850]) / (12400 / 2),
        rs=30,
        output='ba'
            )

b, a
Out[10]:
(array([0.04008421, 0.04822953, 0.04776774, 0.04822953, 0.04008421]),
 array([1.        , 1.66028928, 2.34267679, 1.39001415, 0.70301997]))

Seems like the filter order is twice larger than expected. A similar task using Matlab's "designfilt" performs with the expected number of coefficients. What I am missing? is the filter order a fluid definition?
For lowpass, the number of the coefficients is as expected...

There is a bit of terminology confusion here.
Strictly speaking the order of the filter is the number of poles, so in your case you designed a 4th order filter.
Many software packages use the convention for bandpass/stop filters to specify the order of each slope to make it consistent with lowpass/highpass design calls. That results in a filter of twice the order.
Interestingly enough, even Matlab is inconsistent about that. designfilt('bandpassiir','FilterOrder',2, ...) will give a second order filter. butter(2,[w1 w2],...) will give you a 4th order filter.


Given local responses by a bank of equally spaced (log-)Gabor filters, how can we estimate the response of an intermediate-scale filter? Consider a grayscale image convolved with a bank of 2D wavelet quadrature pairs – in my case, log-Gabor filters. I have eight filters. For simplicity, let's say they are all vertically oriented, and an octave apart in frequency.

At some image location $(x, y)$ I now get 16 response values (one for each even and odd filter at eight different frequencies):

I would now like to estimate the response at an intermediate frequency, say, at three-and-a-half octaves. Can that be done?

Intuitively, I would expect the local distribution of even and odd frequencies to be something to which I can (within obvious sampling limits) fit a function, allowing me to estimate intermediate responses. Is this correct? Can such an expression be derived from the filter kernel itself?

I have since found this paper, which provides a nice framework for how to think about this problem.


Designing FIR filter I'm asked to design FIR filter, an ideal low pass filter.

I'm given with pass band frequency (wp) and stop band frequency (ws). Now, I have to choose Ns points from stop band and Np points from pass band. Taking the inverse will give me impulse response.

The problem is, Ns points are equally spaced and Np points are equally spaced. But the spacing of the points in pass band and stop band differ. How do I find the filter coefficients?

Thank you.

It is not a good idea to have Np and Ns as constraints, for the exact reason you mention- it will usually cause unequal spacing.  Instead, make N, the total number of points, a constraint, and then you will have the equal spacing that you need.


Digital to analog aliasing or mirror query... DAC can output negative frequency? Using a clock Fs = 1 MHz.

The Digital to Analog (DAC) can make frequencies upto 500 KHz.

500 kHz to 1 MHz is an alias? Or is it called a mirror?? Is aliasing a concept only on the ADC and sampling process.

Does -500 to 0 exist as an alias or mirror too... But surely a DAC can't make negative frequencies? Texas Instruments tutorial on DAC has the attached screen shot showing a DAC outputting a negative frequency.

Is anything different for a DAC with complex samples such as on an ADC complex sampling allows frequencies upto Fs instead of Fs/2...!

Block diagram in this link

I think you are confused by negative frequencies and what they mean so let me add this explanation.

When you see a spectrum that contains "positive" and "negative" frequencies, each of the frequencies are of the form:

$$e^{j\omega t}$$

Where $\omega$ is the frequency (in this case angular frequency as $2\pi f$ with f being the frequency in Hz.

The general form of $e^{j\phi}$ is a phasor on the complex plane with magnitude one and angle $\phi$, so $e^{j\omega t}$ is simply a phasor with unit magnitude rotating with rate $\omega$. If $\omega$ is positive, the rotation is counter-clockwise: a positive frequency, and if the rotation is negative, the rotation is clockwise and represents a negative frequency. So the concept of positive and negative frequencies only applies to frequencies of the form of the complex exponential and is the direction of rotation.

Consider Euler's identity for a cosine that demonstrates the concept of negative frequencies well:

$$cos(\omega t) = \frac{1}{2}e^{j\omega t} + \frac{1}{2}e^{-j\omega t}$$

Here we get two phasors on the complex plane, one rotating clockwise, one counter-clockwise, and adding the two results in a signal that always stays on the real axis going back and forth sinusoidally.

So therefore a simple cosine wave with a constant frequency would contain one positive (exponential) frequency and one negative (exponential) frequency. We also see from this that ANY real signal must have positive and negative frequencies, AND the negative frequencies MUST be of the opposite phase. Like the cosine example above, each phasor component spinning in a counter-clockwise direction (the positive frequencies) MUST have an equal magnitude phasor spinning in the opposite direction (negative frequency component) in order for the sum to stay on the real axis (which would then be a real signal with no imaginary components. This is said to be "complex conjugate symmetric" when the negative frequency components have equal magnitude and opposite phase compared to the positive frequency components.

In order to IMPLEMENT a signal that does not have a spectrum with this property (positive and negative frequencies not complex conjugate symmetric), you would need to have TWO values in implementation, one to represent the real number and one to represent the imaginary number (or you could do magnitude and phase, either way you need two signals). You see this in the baseband part of your diagram where the real is represented as "I" for "in-phase" and the imaginary is represented as "Q" for quadrature. So in that part of your diagram we are indeed representing a complex signal, and the spectrum at that point can be asymmetric.

I notated your diagram to demonstrate this further. I will call the IQ signal the "baseband signal" and because it is indeed complex, it's spectrum can be assymetric as I have drawn. This complex baseband signal is then digitally upconverted (frequency translated) by multiplying it by $e^{-j\omega t}$ since  $e^{-j\omega t} = cos(\omega t)-jsin(\omega t)$ and then taking the real part:

$$(I_1+jQ_1)(I_2-jQ_2)$$
$$= (I_1I_2+Q_1Q_2)+j(I_1Q_2-I_2Q_1)$$

Taking the real part of this results in the form in your block diagram prior to the DAC:

$$(I_1I_2+Q_1Q_2)$$

Given as:

$$I cos(2\pi f_{IF}t) + Q sin(2\pi f_{IF}t)$$

With $I_1$ and $Q_1$ as the I and Q inputs to the digital multiplier and $I_2$ is the cosine waveform into the upper digital multiplier and $Q_2$ is the sine waveform input into the lower digital multiplier in the diagram)

So the signal where we maintain two datapaths (truly a complex signal) will have an asymmetric spectrum while the signal where this is only one datapath must be real and therefore must have a spectrum that is complex conjugate symmetric. I like this diagram as shown as it clarifies that the first signal I circled is indeed complex and has an asymmetric spectrum, while the second one I circled, while it looks like the spectrum is asymmetric if you just considered the positive frequencies alone, is not as we are referring to symmetry between the positive and negative frequency axis.

So we see how negative frequencies exist throughout the entire block diagram, and it is not that the DAC is creating them.  The baseband I and Q signal coping out of the encoder is complex and has an assymetric spectrum, this gets digitally upconverted (and is still complex until summed), the output of the summer is now one signal path so must be real and therefore now has a symmetric (complex conjugate) spectrum.  The DAC will convert this to the analog domain and the bandpass filter will filter out the images of the DAC (positive and negative--it is a real filter so it too is complex conjugate symmetric in its passband). So assuming the BPF is centered on $f_{IF}$, the analog spectrum would look identical to the spectrum prior to the DAC, extend to $\pm \infty$ while the unique DAC input spectrum extends only to $\pm f_s/2$.

The "image" shown  in the diagram are NOT the images due to digital sampling, but the image due to the analog upconversion. This is a real multiplier and as such its output will contain the sum and the difference of the input (review the expansion of $cos(\alpha)cos(\beta)$ to see this), so we get the terms RF+IF and RF-IF at both postiive and negative frequency locations. (However you can also see that we are just frequency translating the same spectrum at the output of the BPF by convolving it with the two frequencies of the RF sinusoid.)

There is more not detailed about the DAC upconversion that I didn't cover but you should be aware of including the option to pass a higher Nyquist zone by setting the bandpass filter on one of the DAC images, and that the DAC itself will introduce a Sinc roll--off in frequency due to the staircase reconstruction.


Estimate Filter Coefficients from the Result of Linear Convolution with a Known Signal If I have samples of input say x(1:500) and it passes through FIR filter with 9 taps and some unknown coefficients. The output y(1:508) is also known. The aim is to estimate the filter response in frequency domain. I know using Convolutional theorem $\hat{H} = (X^H*X)^{-1}(X^H*Y)$, but this requires fFT of full 500 samples, However I can only do fft of 128 samples at a time. How to estimate the filter response in frequency domain ?

The idea here is to build the problem in its Matrix Form.
We have the filter $ h $ which is represented by the matrix $ H $ to represent Linear Convolution operation:

$$ H = \begin{bmatrix}
{h}_{1} & 0 & 0 & \ldots &  & 0 \\
{h}_{2} & {h}_{1} & 0 & \ldots &  & 0 \\
\vdots &  & \ddots &  &  & \vdots \\
0 & \ldots & {h}_{m} & {h}_{m - 1} & \ldots & {h}_{1}
\end{bmatrix} $$

This means that the Linear Convolution is given by:

$$ \boldsymbol{y} = H \boldsymbol{x} $$

This is the equivalent of vY = conv(vX, vH, 'full') in MATLAB.
It is clear that for $ \boldsymbol{x} \in \mathbb{R}^{n} $ and $ \boldsymbol{h} \in \mathbb{R}^{m} $ (Namely the filter has $ m $ taps) then $ H \in \mathbb{R}^{ \left( n + m - 1 \right) \times n } $.

The above matrix $ H $ is Toeplitx Matrix.
You may be able to generate it using GenerateToeplitzConvMatrix() I added to my answer repository.

Bu we're after the coefficients of $ \boldsymbol{h} $. The nice thing is we can write the above as:

$$ \boldsymbol{y} = X \boldsymbol{h} $$

Where $ X \in \mathbb{R}^{\left( n + m - 1 \right) \times m} $ is still a Toeplitz Matrix since the convolution (When shape is full in MATLAB notation) is commutative.

So once we have the matrix $ X $ we can solve the above equation like any other linear system:

$$ \hat{ \boldsymbol{h} } = \arg \min_{\boldsymbol{h}} {\left\| X \boldsymbol{h} - \boldsymbol{y} \right\|}_{2}^{2} = {\left( {X}^{T} X \right)}^{-1} {X}^{T} \boldsymbol{y} $$

But since $ X $ has a special structure, being Toeplitz, there must be a way to exploit this, right?
Indeed there is, solving it using the Levinson Recursion. You'll be able to find a solver in the code repository of this question as well (Also see my Levinson Recursion Implementation). The idea is similar to Hilmar's solution.

But I had another idea. Toeplitz Matrix are a generalization of Circulant Matrices. If $ X $ was a Circulant Matrix we could use the highly efficient fft() function to solve the problem.

So the question is how to implement Circular Convolution using Linear Convolution.
Namely, we want to adjust the matrix $ X $ and the output $ \boldsymbol{y} $ so the will represent Circular Convolution between $ \boldsymbol{x} $ to $ \boldsymbol{h} $.

Let's say we know how build those (Let's call them $ \bar{X} $ and $ \bar{ \boldsymbol{y} } $) then we have:

$$ \bar{ \boldsymbol{y} } = \bar{X} \boldsymbol{h} $$

Yet since $ \bar{X} $ it diagonalizable by the DFT Matrix then (See Eigen Values of Circulant Matrix):

$$ \bar{X} = {F}^{H} \Omega F \Rightarrow \bar{ \boldsymbol{y} } = {F}^{H} \Omega F \boldsymbol{h} \Rightarrow \boldsymbol{h} = {F}^{H} {\Omega}^{-1} F \bar{\boldsymbol{y}} $$

Where $ F $ is the DFT Matrix and $ \Omega $ is a Diagonal Matrix of the Eigen Values of the matrix $ \bar{X} $. Those Eigen Values are basically the DFT Transform of $ \boldsymbol{x} $.

So basically we can solve this using the DFT (fft()). In MATLAB code it will be something in the lines of:

vH = ifft(fft(vYBar) ./ fft(vX), 'symmetric');

The full MATLAB code is available on my StackExchange Signal Processing Q45879 GitHub Repository.
It includes implementation of all 3 methods and one can see that each yields a perfect solution as long the system is well defined.

Resources

Circulant Preconditioned Toeplitz Least Squares Iterations.
Conjugate Gradient Methods for Toeplitz Systems.
Scientific Applications of Iterative Toeplitz Solvers.
Iterative Methods for Linear Systems with Matrix Structure.
Least Squares and Total Least Squares Methods in Image Restoration.
Fast Iterative Methods for Least Squares Estimations.
A Weakly Stable Algorithm for General Toeplitz Systems.
Algorithms for Toeplitz Matrices with Applications to Image Deblurring.
A Superfast Toeplitz Solver with Improved Numerical Stability.
Super Fast Solution of Real Positive Definite Toeplitz Systems.
Fast Solution of Toeplitz and Cauchy Like Least Squares Problems (Code in Giuseppe Rodriguez Page - TPLS).
A Super Fast and Stable Toeplitz Direct Solver via Randomization and Hierarchically Semi Separable (HSS) Methods (Code in Jianlin Xia Page - HSS Software Package).


Fourier series approximation: DC component and fundamental frequency In the linked image below, what is meant by plotting the DC component and fundamental frequency for a Fourier series approximation?

For dot point 1 does it want me to graph just the DC component basically or for the first harmonic? Basically what is meant by fundamental frequency?

Plotting the DC component and only the fundamental frequency means, you should plot the harmonic approximation that uses only $a_0, a_1, b_1$ (i.e. the DC component plus sines and cosines of the first frequency (which is the fundamental frequency). For the second point, you should plot the reconstruction with $a_0, a_1, b_1, a_2, b_2$.

Exactly this problem is treated, including code, at my site.


Why is the Fourier transform of a Dirac comb a Dirac comb? This doesn't make sense to me, because the Heisenberg inequality states that $\Delta t\Delta \omega$ ~ 1.

Therefore when you have something perfectly localized in time, you get something completely distributed in frequency. Hence the basic relationship $\mathfrak{F}\{\delta(t)\} = 1$ where $\mathfrak{F}$ is the Fourier transform operator.

But for the Dirac comb, applying the Fourier transform, you receive another Dirac comb. Intuitively, you should also get another line.

Why does this intuition fail?

I believe that the fallacy is to believe that a Dirac comb is localized in time. It isn't because it is a periodic function and as such it can only have frequency components at multiples of its fundamental frequency, i.e. at discrete frequency points. It can't have a continuous spectrum, otherwise it wouldn't be periodic in time. Just like any other periodic function, a Dirac comb can be represented by a Fourier series, i.e. as an infinite sum of complex exponentials. Each complex exponential corresponds to a Dirac impulse in the frequency domain at a different frequency. Summing these Dirac impulses gives a Dirac comb in the frequency domain.


Improvement of SNR by sampling different number of cycles using FFT of a fixed length singal My question relates to Fourier transform of signals with different lengths in time but the same fixed length of the digitized signal with different sampling rates.

Specifically, how does the number of cycles recorded effect the SNR of the FFT spectrum of a signal given a fixed number of NFFT points(the number of samples recorded is kept constant)?

My first signal(square wave of frequency 500 Hz) was sampled at a sampling rate of 100 kilosamples/sec and the signal was sampled for 10 sec, so a total of 1000001 samples.

The second signal(square signal of frequency 500 Hz) i used was sampled for 0.02 seconds at a sampling rate of 50,000 kilo-samples/second. So again the data vector length was 1000001 samples.

The FFT spectrum that i got for signal 1 was much sharper and had a high SNR  compared to signal 2(second signal). intuitively it makes more sense that the signal with more cycles would give a better result.

Is this correct ? more importantly can anyone help me in finding a quantitative relationship for this behavior?

your help will be much appreciated.

Thanks.

Ok, so, what you need to understand is that FFT has particularities that if not taken into account produces erroneous conclusions out of its results. Firstly, both of your signals have the same SNR. Why? The way you described the problem led me understand that you tried to put two identical signals into different process paths to understand why they came out differently. Am I right? So, what are you really trying to ask is “Why both paths led to different results?”. If that’s the truth, allow me:
- FFT can generate lobes around the real frequency in many cases: 1. when you are trying to see a frequency that is not part of the indexes when you selected the FFT lenght; 2. When you don’t use a proper window function prior to FFT analysis; 3. When you have a small window relative to the frequency you’re trying to see and you produce a discontinuity. This is because the FFT algorithm considers that everything you have in your analysis window is a periodic pattern. Probably the 3rd reason is letting you believe that your second measurement has a higher SNR when the truth is your window has only 10 periods of your frequency, and it’s much likely that you have a discontinuity between the beginning of it and the end. If you fix this discontinuity you should see a clear FFT in your second case.


Noise problem when transforming and quantizing difference between images I have a program which does forward and inverse DCT and quantization and dequantization using JPEG matrix. If I apply it to a single image it works fine, but if I try to apply it to a difference between current and previous frame it gets extremely noisy. Album with examples ("high frequency components" line in descriptions means any component starting from anti-diagonal and beyond it).

What is the correct way to compress the difference between two images?

uncompressed

Compressed, high frequency components removed

Difference between this frame and previous one was coded. High frequency components kept

Difference between this frame and previous one was coded. High frequency components removed

JPEG was developed for photographic content, not for difference-of-pictures content. Therefore, the fact that it's not overly applicable to the difference between successive frames isn't surprising.

In your specific quantization approach, where you eradicate high frequencies, that's actually even pretty intuitive: The difference between successive frames are probably small, and fine-grained. Of course, when you erase high frequency content, you lose most of the difference information and get a large error. So, I'd call this an indication that your high-frequency deletion works, and little more.

So, wrong approach to video coding; you could use JPEG for I-frames and code the differences to previous frames with some sort of linear predictor. Or, look at existing video codecs, and adapt some of the methods used there - in the case of grass, which presumably is blowing in the wind in a pretty uniform direction with smaller local variations, motion prediction coding might be of special interest to you.


Simulation of wireless communication system In monte carlo simulation of the wireless communications system why we multiply noise by SNRlin( Signal to Noise Ratio in linear). What is the logic behind doing this. Any help in this basic doubt is highly appreciated.

This is done to change the noise power, keeping signal power same, thus effective changing the SNR. If you keep the signal to unit power and modify the variance of noise in accordance with $\frac{1}{\sqrt{SNR}}$. You would effectively get the desired SNR for the  received signal. In Monte Carlo simulations this is done many times for a fixed SNR to get the average bit error rate for that SNR point. Usually there would be a loop over SNR as well to simulate for different SNR points.

You could acheive the same thing by keeping the noise to be unit varainace and multiple by $\sqrt{SNR}$ to the signal part.

We use the linear scale because that is how the SNR actually impacts reception. The receiver signal model is linear in signal and noise.

Sample code:

  $$y = signal + \frac{randn(1,1)}{\sqrt{SNR}}$$

  Or

  $$y = \sqrt{SNR}*signal + randn(1,1)$$

Both have signal to noise ratio as SNR


How to get discrete values from a nosiy data set that have a trend? I have a sensor that produces values over time (appr. 50/sec). I visualized these values (X axis is the time), this can be seen on the picture, measurement A and B - these are two independent measurements.

I'm a total beginner, I don't even know how to phrase the question. I just need to get... I don't know... "discrete" (?) values from this noisy data set.

I need something like this (this is happening say under 10 seconds):

Up... Up... Up... Up - I'm only interested in the "ups".

I was thinking of summarizing the values and when it's positive that's an up and when it's negative that's a down. But I don't know where to start and until when do I have to summarize the values...

or if is there (I'm sure there is) any other methodology that is used for this kind of problems.

UPDATE:

After talking with Dan Boschen and applying what he recommended (moving average), I got lot smoother result:

Still got difficulties to ignore the false movements, but this is a one-step-forward. I'll continue this topic in a different question.

From the comments you are most interested in counting the number of clear up areas over a certain duration of time, and you point out that you will have occasional erroneous negative spikes even during the up intervals.

You obviously want to filter but you also suggested little knowledge in this space so I will keep it very simple in what should provide a reasonable result.

My suggestion is to do a moving average over the minimum duration you ever expect for an up interval. Pass your waveform through the moving average, and then on the filtered output set a threshold detector with hysteresis (two thresholds) such that you increment your Up counter when passing the upper threshold and then allow your counter to trigger again once your lower threshold is set. Adjust your threshold levels for best results based on your actual signal cases.


Generating Music Tones, How to handle when tones overlap? (16bit samples) If i add the two tones then the result if two peaks are added together will be greater then the the sample max (2^15).

For the life of me I couldn't think of the right words to google. Think i have google block.

If you can't tolerate any distortion, then divide by 2 the amplitude of the signals you want to mix.

Please check this question


baffled by fft phase spectrum! A very simple MATLAB experiment:

f = 200;
fs = 1000;
t = 0: 1/fs : 1;
x = cos(2*pi*f*t);
plot(angle(fftshift(fft(x))));

And here's the output:

Now, made a minor change to the above code snippet; reducing the time duration by just 1 sample, as follows:

f = 200;
fs = 1000;
t = 0: 1/fs : 1 - 1/fs;
x = cos(2*pi*f*t);
plot(angle(fftshift(fft(x))));

And the phase spectrum goes totally crazy:

Questions:

In the first plot, I was hoping to see a zero phase at bin 700 which corresponds to the positive frequency of 200 in this example. That doesn't seem to be the case.
Secondly, I don't understand the linear parts of the graph in plot 1. I do appreciate phase components that could potentially exist due to the so-called numerical noise, but then how could that noise be so 'linear' in phase?
In the second plot, why removing just one sample would have such a drastic impact on the phase plot?
Am I doing something fundamentally wrong here?

You aren't doing anything wrong, but you also aren't thinking carefully about what you should expect to see, which is why you're surprised at the result. For question 1, your conjecture is close, but you actually have things backwards; it's numerical noise that's plaguing your second one, not your first one.

Pictures may help. Here's plots of the magnitude and phase for the first trial:

x = Cos[2.0 \[Pi] 200 Range[0, 1, 1/1000]];
fx = Fourier[x];
ListLinePlot[Abs[fx], PlotRange -> All]
ListLinePlot[Arg[fx], PlotRange -> All]

And the second one:

x = Cos[2.0 \[Pi] 200 Range[0, 1 - 1/1000, 1/1000]];
fx = Fourier[x];
ListLinePlot[Abs[fx], PlotRange -> All]
ListLinePlot[Arg[fx], PlotRange -> All]

So what's going on here? The second one is easiest to explain. First, the FFT of the second one has zero amplitude everywhere except at the two peaks visible in the magnitude spectrum; this follows since the definition of the FFT using 1,000 datapoints returns frequencies of the form $k/1000$ for $0\leq k\leq 999$, and hence your signal falls exactly on a frequency bin. As a result, at the other 998 points your signal is entirely due to machine noise due to floating point errors, and hence your phase spectrum is nonsense, because it's literally the phase of pseudorandom numbers.

For the first one, however, the definition of the FFT encompasses frequencies of the form $k/1001$ for $0\leq k \leq 1000$, whereas your signal frequency is $200/1000$, which is not of the form $k/1001$. As a result, your signal is broadened by spectral leakage, and will be nonzero almost everywhere. I will not comment on the physical form of the phase plot, but I will say that it admits a closed analytic form.

In general, I think graphs of phase angle alone are a really bad idea for conveying information, precisely for this reason; first, you can't tell if you're looking at the phase of low-amplitude garbage or actual signal, and second, it's not translation-invariant, and it's easy to get utterly bewildering graphs for simple inputs. Much better, if you're still looking for something which conveys phase information, is a graph which simultaneously portrays phase and amplitude information in the same visual manner, such as a plot where phase is encoded as hue and magnitude is encoded as brightness.

ADDENDUM: Here's a couple pictures from Mathematica which illustrate the principle I stated in the preceding paragraph:

hue = Compile[{{z, _Complex}}, {(1.0 Arg[-z] + \[Pi])/(2 \[Pi]),
Exp[1 - Max[Abs[z], 1]], Min[Abs[z], 1]},
CompilationTarget -> "C", RuntimeAttributes -> {Listable}];
L = 500;
data = Table[Boole[x <= 11 && y <= 11], {x, L}, {y, L}];
Image[hue@
RotateRight[
10 RotateRight[Fourier[RotateRight[data, {-5, -5}]], {L/2, L/2}]],
ColorSpace -> Hue, Magnification -> 1]
Image[hue@
RotateRight[
10 RotateRight[Fourier[RotateRight[data, {-4, -4}]], {L/2, L/2}]],
ColorSpace -> Hue, Magnification -> 1]
Image[hue@
RotateRight[
10 RotateRight[Fourier[RotateRight[data, {0, 0}]], {L/2, L/2}]],
ColorSpace -> Hue, Magnification -> 1]

All three images are 2D Fourier transforms of the same input signal (an $11\times 11$ square of 1's padded with zeros to a length of $500\times 500$), but the inputs have been cyclically rotated by 5, 4, and 0, and 200 datapoints. The magnitude spectra (encoded by pixel brightness) are identical, but the phase spectra are completely different! Phase encoding is done so that 1 maps to red, $i$ maps to green, $-1$ maps to cyan, and $-i$ maps to purple. This is what I mean when I say that phase spectra are non-shift invariant, and thus not amenable to human visual understanding. For example, with a cyclic shift of 200 datapoints, it's completely impossible to tell what's going on in the phase, since it just looks like static, but the input signal is no more complicated than the other input cases.


Why is the backprojection of the Middlesbury depth map "cone shaped"? I wrote some come to backproject depth maps. Testing it with the Middlebury cone set, I get:

First of all it has noise, which I shouldn't get from the ground truth, but more importantly it seems like the points are going further from each other with the depth.

Can anybody explain this phenomenon?

My code:

  float baseline = 0.3;
  float focalLength = 303.0;
  float numerator = baseline * focalLength;

  for (int i=10; i<groundTruthDisparity.rows-10; ++i) {
    for (int j=10; j<groundTruthDisparity.cols-10; ++j) {
      float distance = groundTruthDisparity.at<uchar>(i,j) / 4.0;

      if (distance != 0 && !isnan(distance)) {
        depthMap.at<float>(i,j) = numerator / distance;
      }
    }
  }

  for(int y=1; y<height-1; y++) {
    for(int x=1; x<width-1; x++) {
      float depth = depthMap.at<float>(y,x);

      // backproject to world
      cv::Matx21d homogeneous = cv::normalize((cv::Matx21d(x,y) - principalPoint) * (1.0 / focalDistance));
      cv::Matx31d transformedPoint = cv::Matx21d(homogeneous(0), homogeneous(1), 1) * depth;

      cloudPoints.push_back(transformedPoint);
    }
  }

The depth map looks like this (colored in HSV):

The reason that pixels get wider with depth is that the depth map works in projection space, so each pixel is basically shaped like a frustum, not like a box.  Not sure where the noise is coming from but the depth map doesn't look perfectly clean either.


Autocorrelation of a signal comprised of only the same number In MATLAB I generated a signal and its autocorrealtion function using:

  signal = ones(1,500).*5; % a vector of the number 5

  acf = autocorr(signal,numel(signal)-1));

To my suprise the resulting ACF vector is all NaN.

But If i generate the signal with numbers between zero and one (0.3 in this case),

  signal = ones(1,500).*0.3; % a vector of the number 0.3

  acf = autocorr(signal,numel(signal)-1));

I get the following plot,

I am just confused as to what is going on here. I would think in either case the autocorrelation of the same number will produce just all '1's indicating no decay what so ever? So why do I get NaN when the signal is above 1 but a linear decaying signal when its between 0 and 1.

Thanks

What you see are mainly numerical effects, even though the second result is actually correct, but not according to what Matlab intends to compute. But let's go slow ...

Matlab computes the autocorrelation function by first subtracting the mean of the signal and then correlating it with itself. Then the autocorrelation function is normalized such that its value for lag 0 equals 1 (this is where you get the NaNs from).

What happens in your case is that the signal equals its mean (because it is constant), so by subtracting the mean you generate a zero signal, the autocorrelation of which is zero. Then in the normalization step you divide the zero signal by its first value (which is zero), and you get NaNs. This shouldn't happen in a well written function, so you could complain with Mathworks about it.

What is even worse, is that this behavior is not consistent. You found some scaling values for which you actually get a reasonably looking autocorrelation function. This is really weird but I think I understand what's going on. Try the following. Look at the values of

signal - mean(signal)

This is what Matlab bases its computiation on. What you'll see is that for the scaling values which give you NaNs, this difference is exactly zero (which is the correct result). For the other scaling values (which give a well-behaved autocorrelation function), this difference will be some constant in the order of $10^{-17}$. So you're again left with a constant signal, albeit it much smaller than before. The autocorrelation function of such a constant signal (without subtracting the mean, because this has already been done) is triangular, simply because it is the convolution of two rectangular functions. The normalization gets rid of the $10^{-17}$ factor and you actually get a nice looking autocorrelation of your signal. This one is even correct if you were to compute the autocorrelation without subtracting the mean of the signal (but this is not what the Matlab function intends to do!).

Conclusion: the first result (NaNs) is useless but at least consistent with the definition of the autocorrelation function used by Matlab (even though a division by zero should be prevented in a good program). The second result is not consistent with Matlab's definition of the autocorrelation, but it is - due to numerical errors (!) - correct when the autocorrelation is defined without first subtracting the mean of the signal. Note that both definitions are used in practice so it's always important to check what is actually meant by 'autocorrelation'.


Going from output of D/C to input of D/C $$y_c[n] {\longrightarrow} \boxed{\text{D/C}} {\longrightarrow} y(t)$$
If the following equation describes the IDEAL D/C converter:
$$y_c(t) =  \sum \limits_{n=-\infty}^{\infty} y_d[n]~ \text{sinc}\left(\frac{\pi (t - nT)}{T}\right)\tag{1}$$
$$Y_c(\Omega) = \begin{cases} T~ Y_d(\omega) \Big|_{\omega=\Omega T} & |\Omega| < \left(\Omega_r = \frac{\pi}{T_s}\right) \\0 & \text{otherwise}\end{cases}$$
How to prove that we can find $y_c[n]$ from $y(t)$ using this formula:
$$y_d[n] = y_c(nT)\tag{2}$$
It doesn't seem entirely obvious to me by rearrange equation $(1)$, that it equals $(2)$.

conceptually, if function $f$ represents a D/C converter, and function $g$ represents a C/D converter, and $f$ and $g$ are inverse functions of each other. then:
$$y(t) = f\Big( y_n[n]\Big)$$
$$g\Big( y(t) \Big) = g\bigg(f\Big(y_n[n]\Big)\bigg)$$
$$g\Big( y(t) \Big) = y_n[n]$$

You won't get the discrete $y_d[m]$ by re-arranging the equations, but by sampling the $y_c(t)$ at integral multiples of $T$. When you sample the RHS in equation1 exactly at $t=mT$ instances, you get $y_d[m]$. You need to remember that $sinc(m-n) = 1$ for $m=n$ and $0$ otherwise, and (BTW in equation 1, in $sinc$ argument $\pi$ won't be there). Anyway, take RHS and evaluate the expression at $t = mT$ where $m$ are integers.
$$y_c(t) \big|_{t=mT} = \sum^{\infty}_{n=-\infty}y_d[n] \, \mbox{sinc} \left( \frac{t-nT}{T} \right) \Big|_{t=mT} = y_d[n]$$


System response: LTI system for $x[n] = \sin(\frac{\pi n}{4})$ Consider the LTI system with frequency response
$$H(e^{j\omega}) = \frac{1-e^{-j2\omega}}{1+\frac{1}{2}e^{-j4 \omega}}, -\pi < \omega < \pi$$
Determine the output $y[n]$ for all $n$ if the input $x[n]$ for all $n$ is
$$x[n] = \sin \left(\frac{\pi n}{4}\right)$$

My attempt: $x[n]$ is an eigenfunction of the LTI system, so the output have the form
$$y[n]=|H(e^{j \omega})|\sin\left(\frac{\pi n}{4} + arg(H(e^{j \omega})\right)$$
But I dont know how to determine the phase and modulus of the frequency response with this form. For example, I think that
$$| 1 - e^{-j 2 \omega} | = \sqrt{(1-\cos(2 \omega))^2 + \sin^2(2 \omega)}$$
Analogous, doing for the denominator, I could not simplificate the result.

The answer:

$$y[n]=2 \sqrt{2} \sin \left( \frac{\pi(n+1)}{4} \right) $$

You're definitely on the right track. The way you're trying to solve the problem is the best and simplest. You just need to realize that you need to evaluate the magnitude and phase of the frequency response just for one frequency, namely the frequency of the sinusoidal input signal:

$$y[n]=\left|H(e^{j\omega_0})\right|\sin\left(n\omega_0+\phi(\omega_0)\right)$$

with $\omega_0=\pi/4$ and $H(e^{j\omega})=\left|H(e^{j\omega})\right|e^{j\phi(\omega)}$.

Evaluating the given frequency response at $\omega_0=\pi/4$ gives

$$\frac{1-e^{-j\pi/2}}{1+\frac12 e^{-j\pi}}=\frac{1+j}{1-\frac12}=2(1+j)=2\sqrt{2}e^{j\pi/4}$$

from which the correct result follows.

As a final note, you incorrectly called the sinusoidal input signal an eigenfunction of the LTI system. This is not true, it is the complex exponential $e^{jn\omega_0}$ that is an eigenfunction, because it appears at the output unchanged other than scaled by a complex constant (the eigenvalue, which happens to be the value of the complex frequency response evaluated at the given frequency). This is not true for a sinusoid, which is not simply scaled by the LTI system.


Relationship of Fourier Transform and Spatial Directivity Pattern I am reading over and over the Section 1.2 of most  understandable mic array tutorial for me  .

A sound wave is dependent on time and space which can be observed in the great link.

But I can't understand why the Fourier transform of the spatially sampled data gives the directivity pattern of the microphone array (which can be seen in the section 1.3.2 of the first document)?

Is there any one who can describe this in more computer sciencist way?

In your second link, the time and space dependence of a simple plane wave is given by

Ψ ( x , t ) = A sin ( ω t ± k x )

Just as a Fourier transform of the equation above can be taken from the time domain $t$ to the frequency domain $\omega$, so also can a Fourier transform be taken from the spatial $x$ domain to the wave vector $k$ domain. For the very simple one dimensional plane wave the wave vector reduces to the wavenumber  $k= \frac{2\pi}{\lambda}$ in that first equation.

In general both the position and wave vector are vector quantities and the dependence of a wave on both is given by

$$
e^{j\vec k\cdot\vec r}
$$

The aperture of a microphone array depends on the position of the microphones (a positional quantity), and so the Fourier transform depends on wave vector.

In Section 1.3.2 of your first link, they are doing exactly that - taking the Fourier transform of the aperture from position to wave vector space, which yields the pattern of the beam in terms of the wave vector.

In the example given in the other answer by @Peter K, think of a plane wave propagating away from the microphone array at an angle $\theta$ .The term inside the exponential, is the dot product of the wave vector and the spacing of the microphone array.


Invariances of FFT-based Image-Registration vs. SIFT-Features TL;DR: I don't understand how invariant FFT-based image-registration techniques are to object alterations (scratches, marks etc.) in comparison to SIFT-features.
I want to build kind of a feature-matching-like system to detect and compare objects from the same production line. The objects have a granular surface structure, therefore matching surface-features already yields good results.
Here are two sample images of objects I would try to identify:

The extracted features should be invariant against translation, rotation, illumination and partial object-alterations, because one use-case could be failure analysis.
Until now I have realized the matching-process by using SIFT or ORB and matching the descriptors of a query-image against a Database of previously extracted descriptors. If I then iterate over the whole database and count the number of good features (using ratio-test, RANSAC or alike) I can find the matching object by searching for the maximum number of good features.
However this process takes a long time due to the high count of descriptors. That's why I looked into techniques that would generate a hash-string or some comparable vector. I found this paper which seems to accomplish a very similar goal to mine. It's based on this paper, which seems to be kind of reference-paper in the field of image registration. There are some implementations available for the latter, however they always solve the problem of rotated/translated template matching (and here).
Finding out how to do feature extraction with these libraries (as in [1]) would be possible (I guess), but I don't understand whether these processes are also invariant to object-alterations! How would the extracted fft-features look, if I alter the object (scratches, marks, dents etc.)? Would I still be able to match them properly?
Another idea would be to extend the existing feature-matching with Bag of Features to speed up the matching process. Would this be more suitable? Is there another, "better" approach that I haven't considered yet?

I don't understand whether these processes are also invariant to object-alterations!

They are not.

How would the extracted fft-features look, if I alter the object (scratches, marks, dents etc.)? Would I still be able to match them properly?

Yes. But within certain limits.
The Fourier-Mellin transform is based on the Discrete Fourier Transform (DFT) whose coefficients depend directly on its input image's values. Any difference in some $I(x,y)$ will be propagated to $\mathcal{F}(I(x,y))$ and subsequent stages.
In the Fourier Mellin transform, this propagation is even longer because the transform's output is the normalised Discrete Fourier Transform (DFT) of the cartesian-to-polar followed by log transformed amplitude DFT. That is, two DFTs, one on the image and one on circular paths around the amplitude spectrum of the image. But the end result is the same. Changes in $I(x,y)$ propagate to the transform's output. Both transforms are inversible by the way.
In the context of your application, this would introduce some "fuzziness" around the centroid of your classes (For example "Part is good", "Part is defective"), which would be expected. If the texture differences (dents, scratches, etc) are in the very high frequencies it would not even matter for techniques like the one described in the Ishiyama paper where the match is based on part of the spectrum (low frequencies, or, the overall shape of the part).

The extracted features should be invariant against translation, rotation, illumination and partial object-alterations, because one use-case could be failure analysis.
...
Another idea would be to extend the existing feature-matching with Bag of Features to speed up the matching process. Would this be more suitable? Is there another, "better" approach that I haven't considered yet?

It is impossible to achieve all of these goals just at the level of the features.
First of all, what you are dealing with here is probably a One class classifier. You know what "good" looks like and you want to allow parts that are within very close proximity to this estimation of good.
The second point to clarify is the nature of the images. How do defects present themselves on the image? What does a deffective part look like? What does a good part look like? What are the processes that cause these defects? Are you interested in a "yes/no" decision from an image of a part or are you interested in localising where the defect is on the image?
To give you an example of the kind of reasoning you might want to use here: Suppose that we want to quantify how flat a surface that is produced by a planing machine is. We can illuminate it with a source at a very shallow angle that would cause any dents or pits or scratches, to appear as dark streaks or spots on the surface image. Then, simply by thresholding the image and measuring how many "dark pixels" appear in it gives us an estimate of how smooth that surface is. Notice here, there was no reason for fancy features, etc. Only knowledge of the process and the output we are after.
Finally, illumination differences in an industrial environment are easily controlled. And if you have a challenging subject like metalic parts that can produce reflections, you definitely want to control the lighting. You can apply some sort of equalisation or overall normalisation of the brightness, but it is impossible to control for ALL of the artifacts that reflections might produce (epsecially from a single image).
Hope this helps.


What is the difference between t=n*Ts sampling vs impulse train sampling? I know that if I sampling with impulse train so I get in the frequency plane
X(f)*h(f)
(when x(f) is my signal, * means convolution and h(f) is fourier transform of impulse train).

what the difference between this way to x(t=n*Ts) sampling way?

If I can get explenation in the "frequency plane" it will be great.

I cannot recall a meaningful difference between the ideal impulse train sampling and the uniform sampling relation indicated by the expression $t_n = n T_s$.

Given a continuous-time bandlimited signal $x_c(t)$, when you sample this with an ideal impulse train $\delta_{T_s}(t) = \sum_k \delta(t - k T_s)$, then the relation between the obtained discrete-time sequence $x[n]$ and the continuous-time signal will be given by :

$$ x[n] = x_c(t_n) = x_c(n T_s)$$

Hence, the latter expression is a consequence of the former operation.


Why are there no poles for finite length signals My professor recently said something along the lines of, "finite length signals can't have poles because poles correspond to recursive difference equations which correspond to infinite length signals."

I can't really wrap my head around this. Can someone explain in a different or more simple way why finite length signals can't have poles? Would much prefer a simple, hand-wavy explanation that I can remember easily rather than a rigorous explanation/proof.

Hand-wavy? A rational fraction $$\frac{N(X)}{D(X)}$$
where $N(X)$ and $D(X)$ are polynomials can be written as some polynomial $P(X)$ plus a finite amount of terms of the shape:

$$\frac{a_i}{X-X_i}\,.$$

Let us say that a $0$ pole ($X_i=0$) is harmless, because it is
 just some $X^{-1}$. So  on the one hand, you have the monomials $X^d$, $d\ge0$, from $P$ and possibly a $X^{-1}$. If you replace the $X$ by a $z^{-n}$, you will still have a finite number of $z$ to the power of some integer, which are just delays.

On the other hand, if you have one non-zero pole:

$$\frac{1}{X-X_j} = \frac{1}{X}\left(1+ \left(\frac{X_j}{X}\right)^1+\left(\frac{X_j}{X}\right)^2+\cdots\right)$$
and replace the $X$ by some $z^{-n}$, you will get an infinite quantity of $z$ to the power of some integer, and thus an infinite impulse response (IIR). The complicated part is to show that several non-zero poles do not "cancel" each other. A lot of people confuse IIR and recursive.

Be cautious though that an apparently "recursive implementation" does not always imply non-zero poles: a FIR system with $z$-transform $1-z^{-2}+z^{-4}-z^{-6}$ can be written as:

$$ \frac{1-z^{-8}}{1+z^{-2}}$$

that is, you can rewrite: $$y ( k ) = x ( k ) - x ( k - 2 ) + x ( k - 4 ) - x ( k - 6 )$$ as $$y ( k ) = x ( k ) - x ( k - 8 ) - y ( k - 2 )\,,$$

which is a recursive difference equation, yet with a FIR.


What is a window function in DSP and why do we need it? When designing digital filters we start off with FIR and IIR. Later on we talk about window functions e.g Kaiser and Hamming windows.

I do not understand why we need so many type of window functions described by rather complex mathematical equations and how they relate to FIR and IIR.

The FFT computes the coefficients for a periodic, ever repeating signal. In practical applications, signal segments are finite and not periodic. Just applying the FFT to a segment of a signal essentially wraps the end of the segment around to the start, generating a jump discontinuity. Such jumps result in an undesirable background artifacts in the FFT amplitudes. To make the wrap-around smooth, one can fade the signal to zero at both segment ends by multiplying it with a window function.

If the signal is essentially smooth, then even a simple linear fade to zero may not enough, since also kinks in the signal lead to some background artifacts in the amplitudes, however much smaller than those resulting from jumps. The choice of the windowing function is thus a balance between background artifact reduction and computational simplicity, and also of the length of the segment.

Additional considerations are required for overlap-add schemes, where the shifted window functions have to add to one.


Computer vision, Coupon code detection, sugested algorithms The telecommunication operator in my country has coupons that used to recharge the balance

I want to develop an Android application, that with the camera can take the scratched number and recognize and automatic recharge.

I think finding the number box, and latter preform an OCR or a pattern match in because are numbers only

I'm new in computer vision I set the environment and test some opencv ups, but I'm losing witch algorithms or sequence of algorithms may use to achieve that I want Exist segmentation, feature detection, blob detection, ... template match but the image maybe rotates scale because the use of the camera

And I'm lose I wonder if exist one algorithm to find the 4 extreme points of a coupon in one image from here I know how to proceed until the OCR part I don't want the code to do it of course

What I want is a simple clue, answer like you may use X algorithm to have Y result

Thanks and sorry for my english

You may try an OCR detector like Teseract with OpenCV:

https://github.com/tesseract-ocr

Train the Teseract algorithm with text samples obtained from the numbers you want to detect and it will work.

EDIT: If you want a similar implementation of what you want you may check this:

http://answers.opencv.org/question/10985/how-can-i-determine-the-location-of-lcdled-display-in-an-image-with-opencv-and-perl/


Autocorrelation of filtered cyclostationary random process Background:

When a wide-sense stationary (WSS) random process $x(t)$ is passed through a linear, time-invariant (LTI) filter with impulse response $h(t)$ to produce $y(t)$, the following relationship holds:
$$
S_y(f) = |H(f)|^2 S_x(f) = H^*(f) H(f) S_x(f),
$$
where $S_x, S_y$ are the power spectral densitities (PSD) of $x$, and $y$, respectively, and $H(f)$ is the Fourier transform of $h(t)$ (ordinary frequency, unitary).  In the time-domain, this relationship can be shown to be:
$$
R_y(\tau) = h^*(-\tau) * h(\tau) * R_x(\tau),
$$
where $R_x,R_y$ are the autocorrelations of $x$, and $y$, respectively.

Question:

If $v(t)$ is a wide-sense cyclostationary process with autocorrelation $R_v(t, \tau)$, is cyclic with period $T_0$, and is passed through the same linear filter to produce $w(t)$:
$$
w(t) = h(t) * v(t),
$$
will $w(t)$ also be cyclostationary with cyclic period $T_0$? (Assuming $h(t)$ is BIBO stable.) If so, is there a "simple" relationship that expresses $R_w(t, \tau)$ in terms of $h$ and $R_v$?

EDIT:

For a wide-sense cyclostationary process $v(t)$, the cyclic autocorrelation $R_v^{n/T_0} (\tau)$ is defined as:
$$
R_v^{n/T_0} (\tau) = \int_{-T_0/2}^{T_0/2} R_v(t, \tau) e^{-j n 2\pi t/T_0} dt
$$

A relationship between $R_v^{n/T_0} (\tau)$ and $R_w^{n/T_0} $ would also be acceptable for this question.

Your first two equations generalize to the spectral correlation and cyclic autocorrelation. The most general results are (12) and (13) on my blog post on input-output relations for $n$th-order cyclic cumulants and $n$th-order cyclic polyspectra: signal processing operations and cyclostationary signals.

You can find the generalization of your first equation in several places, including Professor Gardner's book "Statistical Spectral Analysis", Eq (90) in Chapter 11:

$$S_y^\alpha(f) = H(f+\alpha/2) H^*(f-\alpha/2) S_x^\alpha(f)$$

And yes, the output of the filter is also cyclostationary with the same period provided that the transfer function $H(f)$ is such that $H(f+\alpha/2)H^*(f-\alpha/2)$ is not zero for at least some $f$ for which $S_x^\alpha(f)$ is also not zero. In other words, the filter could remove so much of the energy of the signal that the remainder is no longer cyclostationary.


Histogram of Gradient I was calculating which would be faster in computation, Inverse Fast Fourier Transform or Histogram Of Gradient but I am unable to calculate its running time.

What would be the running time of Histogram Of Gradient if the image has 'p' pixels ?

For the sake of putting some numbers to this question, I implemented a basic histogram of gradient from scikit-image (skimage.feature.hog). Here is the timing data for HOG with default parameters applied to the skimage.data.astronaut image in b&w and rescaled to have the given dimensions:

Image dimensions......(102, 102)
3.49 ms ± 26.3 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Image dimensions......(256, 256)
27.5 ms ± 295 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Image dimensions......(512, 512)
123 ms ± 2.88 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
Image dimensions......(768, 768)
284 ms ± 9.18 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
Image dimensions......(1024, 1024)
509 ms ± 8.62 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
Image dimensions......(2048, 2048)
2 s ± 15.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

And here is timing data for the inverse Fourier transform (np.fft.ifft2) on the Fourier transform of the same set of images. (Note that taking the Fourier transform is not included in the timing results)

Image dimensions......(102, 102)
684 µs ± 32.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
Image dimensions......(256, 256)
2.51 ms ± 26.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Image dimensions......(512, 512)
21 ms ± 519 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Image dimensions......(768, 768)
42.5 ms ± 621 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Image dimensions......(1024, 1024)
85.9 ms ± 1.47 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
Image dimensions......(2048, 2048)
381 ms ± 9.18 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Graph of the timing data


An Impulse in Complex Envelope (i.e. baseband) Notation? In complex envelope notation, a modulated RF signal, say v(t), is represented in complex envelope notation as v(t) = Re{ g(t)e^(jwct)}, where wc is the carrier angular freq, g(t) = x(t) + jy(t), and g(t) is the complex baseband signal.  Since x(t) and y(t) are inputs to an in-phase (cos) and quadrature phase (sine) mixers, how would I represent v(t) in baseband g(t) notation, where v(t) = an impulse?  Since a mixer is technically nonlinear, and therefore NOT an LTI system, isn't this impossible to do in a mathematically correct way?
The issue is we have a MATLAB RF simulator for PHY signals, whereas the simulator is designed to work only with RF signals represented in baseband notation (and fc is metadata), and I want to send impulses through the RF channel to clearly see the fast fading reflections.  My current method is to use impulses for x(t), with y(t) = 0, but not sure if this translates correctly, or if this even makes sense.

The OP commented that he has practical reason he would like to pursue a simple bandpass equivalent of an impulse test. Here are some considerations with that in mind:
An impulse in time represents all frequencies, which will ultimately be band-limited by the channel. We could implement the IQ representation of the band-limited channel that could then be transmitted at a real carrier if the carrier frequency itself far exceeds the bandwidth (which is often the case). In this case, the impulse could be implemented as a truncated Sinc function on I and Q to maximize the channel power, or even on just I if that isn't a concern in such a way that the Fourier Transform of the Sinc (which is flat in frequency) sufficiently exceeds the channel bandwidth, which would then be no different than an impulse itself. So this means the ideal solution would be to transmit Sinc shaped RF pulse envelopes in time such that the entire channel is occupied evenly in frequency, which can be well approximated from a truncated Sinc on I and Q at baseband. This approach would maximize the energy in the band of interest while providing the same result as an impulse.
Prior comments below on more typical channel estimation techniques:
Any waveform that occupies a particular bandwidth can be identically represented at any carrier including a carrier frequency of zero (which is DC) as long as we use complex representation to allow for different positive and negative frequencies in the overall bandwidth. This is often much much simpler for modeling and expression given we don't have to represent each cycle of some higher frequency carrier. A mixer as a pure multiplier is non-linear in that it's output frequencies do not match it's input frequencies but as such it is simply a frequency translation with no other modification to the signal. Actual mixers will have other non-linearities generating other frequencies that we control with good design practice as with any other analog component.
The second point is to consider NOT using an impulse for channel estimation. In practice any impulse that can be created will have so much less energy than other techniques often used, and as such will provide far less SNR in the answer. The two most common approaches for channel estimation are swept tones (such as that done in a network analyzer) or psuedo-random (PRN) sequences which emulate a white noise source and as such occupy the channel. The difference is both of these approaches can provide a sustained signal for as long as needed for desired SNR only limited by how stationary the channel is. The approach for channel estimation with a pseudo-random sequence would be ideal for a fast fading channel, with the duration of the test sequence set to be less than the fade rate of the channel. I detail channel estimation using PRN sequences further in these posts:
This shows an equalizer implementation but by swapping input and output you would instead get channel estimation:
Compensating Loudspeaker frequency response in an audio signal
and
How determine the delay in my signal practically
These approaches will return the impulse response in time, from which you see a single dominant impulse (for a channel with no reflections) or multiple impulses (from the multiple reflections, with their corresponding magnitude and phase).


Digital signal derivative I'm new to signal processing and I need your help: I have an array of 128 elements (call it Window) filled with 128 samples taken from a sensor. I was wondering how to find the derivative of this digital signal. I found this formula on wikipedia:

m = F[n] - F[n-1]

where, for me, F[n] is Window[n]. I'd like to obtain an array M filled with all the derivatives. But there's something I miss: The first step would be:

M[0] = Window[0] - Window[-1]

and I haven't Window[-1]. How can I manage this issue? Thank you in advance.

https://en.wikibooks.org/wiki/Digital_Signal_Processing/Discrete_Operations#Derivatives

In theory, M[0] doesn't exist. Your discrete derivative is only defined on [1,127].

As a practical solution, you can define M[0] as a null element of your signal. For example if you expect values between -1 and 1, take M[0] = 0.

In the case of real-time, there is two options:

Either you save the last value of your buffer so you can take M[0] = last for the next ones, and initialize last as 0, or whatever seems coherent with your expected signal;
Either you start computing the derivative at M[1] and your processing will have a latency of 1 sample.


Generating wavelet using amplitude spectra Equation (4) of a recent paper (An improved peak-frequency-shift method for Q estimation) shows how to generate a wavelet using a zero-phase Ricker wavelet.

$|B_1(f)| = |B_0(f)|\exp(-\pi t f / Q)$

In the above, $|B_1(f)|$ is the amplitude spectra of the wavelet that I want to calculate, $|B_0(f)|$ is the amplitude spectra of the reference (Ricker wavelet), $\pi$ is the ubiquitous constant, $t$ is the time separation between wavelets, $f$ is the frequency, and $Q$ is a number (the Quality factor).

I know how to calculate the zero-phase Ricker wavelet in the time domain as $B_0(t)$, and I know how to calculate $|B_0(f)|$ as the square root of the sum of the real and imaginary parts squared of the frequency-domain reference Ricker wavelet.  The real and imaginary parts are found using the FFT.

However, now that I've calculated $|B_1(f)|$, how do I transform $|B_1(f)|$ back into the
time domain?  I know $|B_1(f)|$, but can I find the real and imaginary parts of $|B_1(f)|$ so that I can use the IFFT to compute $B_1(t)$ in the time domain, or is some other information missing and required?

Is the "without phase changes" words in the recent paper (http://www.seg.org/documents/10161/997244/1218.pdf) a clue that might help with the numerical implementation?

So how did the authors generate the wavelets in Figure 1 of the paper?  Is something other than Equation (4) being used here?

Perhaps there is nothing more happening here than an issue with notation. Looking at another paper (Estimation of quality factors from CMP records), the frequency spectrum of a Ricker wavelet is computed using Equation (1):

$B(f) = \frac{2}{{\sqrt \pi  }}\frac{{{f^2}}}{{f_m^2}}{e^{\frac{{ - {f^2}}}{{f_m^2}}}}$

In the above, $f_m$ is the dominant frequency.  However, there doesn't seem to be any complex part of this expression.  Does this imply that the complex part is comprised of all zeros?  To numerically evaluate this equation, how do I compute the $f$ variable?  For example, using Matlab, $f$ can be computed using the colon operator:

f = fmin:deltaf:fmax

But what is fmin and fmax?  Now $B(f)$ can be computed and taken back into the time domain using the IDFT (or the IFFT), but what if I want a real time-domain discrete signal for the $B(t)$ function?

Equation (5) of the same paper gives an expression that I think can be used to calculate the Q factor wavelet (note that this looks the same as the equation at the very top of this post):

$B(f,t) = B(f){e^{\frac{{ - \pi ft}}{Q}}}$

However, what is meant here by the "amplitude spectrum"?  How do I get the real and complex parts of this spectrum?

A time domain version equation of the Ricker wavelet can be found here: Wikipedia link, but then what is the time domain expression used to calculate the Q factor wavelet?  As an aside, it would have been easier to have the time domain version of:

$B(f,t) = B(f){e^{\frac{{ - \pi ft}}{Q}}}$

But what is the "amplitude spectrum" in this context, and how can I compute the wavelet in the time domain?

The complex result of a DFT of a signal codifies amplitude and phase information. When you take a magnitude of the resulting complex vector, you automatically throw away any phase information.

In so many words, the DFT result is informing you of what amplitudes and phases of all its frequency components are required to reconstruct your signal. If you strip away phase information you are lining up all your complex sinusoids and then adding them up coherently. If you lose phase how can you reconstruct your signal?

That being said, what you can get in this case is the autocorrelation function of $b_1(t)$, by taking the IDFT of $|B_1(f)|^2$, but you need phase information to completely reconstruct the original signal.


Adaptive Thresholding for liver segmentation using Matlab I need to segment liver from an abdominal CT image using Adaptive Thresholding. But I get the entire foreground separated from the background alone. I need only the liver part of the foreground separated.
Check out the pdf file in http://www.ijcaonline.org/casct/number1/SPE34T.pdf
I need an output similar to that shown in Figure 6.

I attach my coding here. Kindly help me out.

%testadaptivethresh.m
clear;close all;
im1=imread('nfliver2.jpg');
bwim1=adaptivethreshold(im1,11,0.03,0);
figure,imshow(im1);
figure,imshow(bwim1);
imwrite(bwim1,'at2.jpg');

function bw=adaptivethreshold(IM,ws,C,tm)
%ADAPTIVETHRESHOLD An adaptive thresholding algorithm that seperates the
%foreground from the background with nonuniform illumination.
%  bw=adaptivethreshold(IM,ws,C) outputs a binary image bw with the local
%   threshold mean-C or median-C to the image IM.
%  ws is the local window size.
%  tm is 0 or 1, a switch between mean and median. tm=0 mean(default); tm=1 median.
%
%  Contributed by Guanglei Xiong (xgl99@mails.tsinghua.edu.cn)
%  at Tsinghua University, Beijing, China.
%
%  For more information, please see
%  http://homepages.inf.ed.ac.uk/rbf/HIPR2/adpthrsh.htm

if (nargin<3)
    error('You must provide the image IM, the window size ws, and C.');
elseif (nargin==3)
    tm=0;
elseif (tm~=0 && tm~=1)
    error('tm must be 0 or 1.');
end

IM=mat2gray(IM);

if tm==0
    mIM=imfilter(IM,fspecial('average',ws),'replicate');
else
    mIM=medfilt2(IM,[ws ws]);
end
sIM=mIM-IM-C;
bw=im2bw(sIM,0);
bw=imcomplement(bw);

My modified code for testadaptivethresh.m

clear;
im=imread('nfliver7.gif');
figure,imshow(im)
bwim1=adaptivethreshold(im,300,-0.15,0);
bw=bwareaopen(bwim1,3000);
se=strel('diamond',4);
er=imerode(bw,se);
bw1=bwareaopen(er,3000);
er1=imerode(bw1,se);
bw2=bwareaopen(er1,1000);
fi=imfill(bw2,'holes');
figure,imshow(fi)

op=uint8(fi);
seg=im.*op;
figure,imshow(seg)
imwrite(seg,'sliver7.jpg');

I saw the link of the paper (of SS Kumar) you referred and the link where you got the code (HIPR) are two different algorithms - even though both sounds like Adaptive Thresholding

First i would like to tell you the difference.

In the HIPR method, a general assumption is essentially 2 class level image - i.e. foreground and background. In any 2 class thresholding, one expects 2 peaks or regions within the image histogram, specifically background vs. forground, text vs. white paper. If you somehow found an optimal valley point in the histogram - you get cleanest division. Here is how the histogram may look like :

However, this valley point may be slightly shifting around the locally. There are good examples of lighting variations given there. Hence, this optimal valley point exists everywhere but slightly varies spatially, hence a universal threshold would fail. Hence, the valley point (threshold) is computed on every local region.

SS Kumar's paper's method and more specifically the class of images you are dealing with, is multi-class (multiple objects each with different intensity band and spread). In this cases, histograms are multi-modal, i.e. it has many peaks and valleys and presumably each peak correspond to different object, yet it may be even more complex.

The histogram might just look like this: (this is the same image as in paper).

In this case, the above 2 class approach will simply fail because there is no one good valley. Which is why your first image you posted looks like black/white dots all around.

The meaning of Adaptive Thresholding here, implies that you need to identify the correct peak and the band of gray scale which covers most intensities of the lever and other objects are in stark contrasts which allows

What should you do?

First of, if it is compulsory to use Adaptive thresholding, find the histogram and see what intensity range and then for a threshold of left or to the right are the intensity boundaries which pixels should be discarded.

Alternatively you can use Region Growing or split and merge algorithm. Refer to this question for some information: What segmentation methods can be used for simple images?


Need information about following channel model The following research paper is for blind channel equalization and the author has modeled the channel as 7 tap weights

https://ieeexplore.ieee.org/document/1096877

On page 3 it has shown the channel tap weights which they have used it. I want to know information about the channel? Is it the multipath fading channel or what kind of channel they are using? If somebody can help? and can I simulate the following channel in MATLAB,if yes then how?

The channel impulse response is given in figure below:

The channel impulse response represent the weighted gains and phases at different delays, so yes this is a multipath channel. You can simulate the channel (in Matlab as you requested, what I show below was done in Octave with compatible results) by first creating the complex channel weights:

channel = [c1, c2, c3, c4, c5, c6, c7];

where $c1 = -.005 - j0.004$, $c2 = 0.009 + j0.030$, etc

The following are examples of inspecting and using the channel:

filter

Use the filter command to determine the complex baseband received signal from the complex baseband transmitted signal (analytic signal):

rx = filter(channel, 7, tx)

Where rx is the received signal and tx is the transmitted signal. The denominator factor of 7 is somewhat arbitrary and is just a scaling.

freqz

Use the freqz command to plot the magnitude and phase frequency response of the channel.

freqz(channel)

The horizontal axis is normalized angular frequency from DC to $2\pi$ where $2\pi$ represents the sampling rate used. (So if the tx data was sampled at 1 MHz for example, the right side of the plot would be 1 MHz). Because it is a complex spectrum the full frequency range from DC to the sampling rate is unique as shown.

grpdelay

Use the grpdelay command to determine the group delay of the channel, which is the delay of the signal in time or in samples at each frequency. In the phase vs frequency plot above it appears the channel may be linear phase, however a group delay plot (group delay is $d\phi/d\omega$) would give us a better indication since it will remove the linear trend.

grpdelay(channel, 1, 512, 'whole');
xlabel('Normalized Frequency (rad/sample)')
ylabel('Delay (samples)')
title('Group Delay')

From the above results we see that the channel distortion is not significant or challenging both in magnitude (+/- a few dB) and group delay (+/- half a sample), in comparison to a highly frequency selective channel that would have deep nulls throughout the passband magnitude response. The delay distortion is more pronounced if a lower sampling rate is used, so depends on how many samples per symbol there are (The delay can be converted to seconds by dividing by the sampling rate).

See this post on implementing an LMS Equalizer to compensate for this channel in the receiver (using the Wiener-Hopf equations): Compensating Loudspeaker frequency response in an audio signal .

For "blind" equalization using this approach, I would use receiver decisions to recreate the expected transmitted signal (including any pulse shaping etc). I did not read your reference link but you could do that as a comparison if that wasn't the technique used directly. A comparison should include performance under lower SNRs.


Sequence expansion by zeros and interpolation - does it insert additional frequencies? I am struggling with understanding the consequences of oversampling on the frequency spectrum of the signal.

If I understand correctly, with an oversampling rate of 8X we insert 7 new values for each value measured. In case we keep this new values at 0 then we are not inserting any new frequencies to the signal, so the frequency spectrum should stay the same.

On the other hand, if we interpolate the new values eg by averaging the adjacent ones, then we MIGHT generate some new frequencies, that might appear in our frequency spectrum. This new inserted freuquencies will be of higher magnitude and should be filtered out if we want to keep the signal as "clean" as possible.

Am I correct with my understanding of oversampling and interpolation? Thanks in advance for any answer!

The key to understanding what inserting zeros does is to understand two things: what samples represent in the time domain (because we want to insert zeros in the time domain), and what they represent in the frequency domain (because we want to know what it did to the spectrum).
First, sampling is a type of modulation (PAM—Pulse Amplitude Modulation), equivalent to multiplying an impulse train by our analog signal, and creates images in the frequency domain. When converted to digital sample values, we have PCM (Pulse Code Modulation), a common term for digital audio. Here is an example spectrum represented by our samples; the original analog signal's spectrum is shown in green, and images in red:

These images are the price we pay to represent the analog signal as samples. As such, the usable bandwidth is from 0 Hz up to (but not including) half the sample rate. Above that is a backwards image of our original spectrum, with images repeated around multiples of the sample rate. It's OK, we remove the images when we convert back to analog, using a DAC's lowpass filter.
In the time domain, samples represent impulses. They are instantaneous values taken at a constant interval. Inserting zeros changes nothing except what we consider the sample rate.
For instance, sample a signal one time per second. That represents an impulse train, so consider playing it back raw, as an impulse train.
Now consider placing a zero-valued sample between each of the original samples. Consider playing it back as an impulse train but at twice the original rate, two times per second.

Can you see that nothing has changed in the signal, except the sample rate? Likewise, if we look at the spectrum, nothing has changed in the frequency domain—this is obvious, since the time domain signal has not changed.

However, our usable bandwidth has doubled. The first inverted image now lies in our usable band, now shown in green. It will not be removed by the DAC when played back through a DAC at the new, higher rate. And it will be a problem with any non-linear processing in the digital domain.
That is why we follow the zero-insertion with a lowpass filter (or combine the two steps for efficiency). Here it is again after proper filtering, below half the original sample rate:

So, the answer is that nothing changes in proper integer sample rate conversion by zero-insertion, and the result is as perfect as the lowpass filter used to clean up the exposed images.


What will the impulse response of a matched filter look like if the input is complex? If $h(t)$ is the impulse response of a filter matched to a signal $s(t)$, I read that $h(t) = ks(t_o - t)$.

But what if the signal is complex? I went through the derivation of the matched filter and saw that just like how the inner product of two vectors is maximum if the vectors are "parallel". I think $H(\omega)$ has to be proportional to $S^*(\omega)e^{-\omega t_0}$ for the SNR to be maximized. This is a generalization of the Cauchy–Schwarz inequality that I don't fully understand.

So from this we can say, $H(\omega) = S^*(\omega)e^{-\omega t_0}$

To get RHS,

frequency shifting,
$$s(t-t_o) \iff e^{-j\omega t_0}S(\omega)$$

taking conjugate on both sides,
$$s^*(t-t_o) \iff e^{j\omega t_0}S^*(-\omega)$$

scaling with $-1$,
$$s^*(t_o-t) \iff e^{-j\omega t_0}S^*(\omega)$$

So can I conclude $h(t) = s^*(t_o-t)$??

To apply this general formula to a specific example, consider this question from the 2013 exam of ISRO,

From what I learned, can I say that the answer is option D?

The impulse response for the given complex signal $s(t)$ is complex conjugate time reversal of $s(t)$, i.e. the real part of $s(t)$ simply gets time reversed (flips across y axis) so A B or D, while the imaginary part is flipped across both the y and x axis. C or D.

$$h(t) = s^*(t_o-t)$$

$$h(t) = (x(t_o-t)+iy(t_o-t))^*$$

$$h(t) = x(t_o-t)-iy(t_o-t)$$

So my question is, what is the impulse response of a matched filter if the input is complex and why? Like please point out any mistakes in my logic, I especially would love to gain some intuition on why $H(\omega)= kS^*(\omega)e^{-\omega t_0}$

The impulse response doesn't change because it is the output of the filter when the input is an impulse.

When the complex input is applied to a complex filter the output is still determined by the convolution of the input and the impulse response of the filter. The addition/integration and multiplication operations that take place inside a convolution operation have to account for complex numbers.

I think it is difficult to gain intuition for the result of a matched filter in the time domain. Instead look at it from the frequency domain.
\begin{align}
 Y(\omega) &=&  H(\omega)S(\omega) \\
 & = &  e^{-j\omega t_0}S^*(\omega) S(\omega)\\
 & = &  e^{-j\omega t_0} ||S(\omega)||^2
\end{align}

In the case where the envelope of $S(\omega)$ is rectangular (as in your example) then $||S(\omega)||^2$ is equal to a rectangular pulse. The output then be a time delayed version of a sinc() pulse. The wider the frequency range of $S(\omega)$ then the narrower the main pulse of the sinc() pulse.

If you are trying to gain insight into how a matched filter is derived - one insight I can provide is that you have to understand that the Signal to Noise ratio (SNR) is defined differently at the input of the filter, than the SNR that we are trying to maximize at the output. At the output we are interested in the maximum value of the output signal, while at the input the SNR is measured using the signal's energy.