Is it possible to use MRC with DS-CDMA-MIMO I'm wondering how to use the Maximum Ratio Combining with MIMO-DS-CDMA, I have used the MIMO with DS-CDMA using Walsh code, then at receiver I couldn't recover the transmit signal using the MRC.  here what what I did:

x = [0.7+0.7i; 0.7-0.7i; 0.7i+0.7i; -0.7+0.7i];    %generate signal x size 4x1
h = (randn(4,4) + 1j.*randn(4,4))/sqrt(2);      %Generate 4x4 channel
code = hadamard(4)                                  %generate walsh code
x_s = kron(x,cod(:,1))                              %spread the signal x with one walsh code
r = reshape(h*reshape(x_s,4,[]),[],1);              % the received signal with AWGN

suppose I have perfect knowledge channel at receiver side, it means h is known.  how can I use MRC for r in order to recover the transmitted signal? Is that possible?

thank you

OK, let me explain that for you. First, if you want to understand anything in engineering, it's recommended to write it in mathematical form, then try to solve it theoretically and build its code accordingly. Second, regarding your question, is it possible to use MRC for MIMO - CDMA,   Yes, that's possible, why not?

Now, The code you provided is not right. Why ? all simply, first   check when using simple MIMO without spreading spectrum. You will see that first element in x, which is 0.7+.7i, is multiplied by the first column of channel H and so on. In case if you are using spreading spectrum, it means that instead of having one sample 0.7 + 0.7i, you will have N samples, where N is the code length. it means you will have N samples in the same duration of symbol before spreading.  So, all samples should be multiplied by the first column of H.

I think it's clear now, right?  so you should just now check how to write the above answer in matlab. . .  for example check this command :

r = [];
for kk = 1 : size(h,1)
                      r_1 = kron(x_s(:,kk),h(:,kk));
                      r = [rec r_1];
                  end

r = sum(r,2);

Good luck.


MATLAB: $\tt fft$ and $\tt ifft$ scaling In MATLAB, the outputs of the fft and/or ifft functions often require additional processing before being considered for analysis.

I have heard many differing opinions on what is correct:

Scaling

Mathworks states that fft and ifft functions are based on the following equations:
\begin{align}
X[k] &= \frac{1}{1} \cdot \sum_{n=1}^{N} x[n] \cdot e^{\frac{-j \cdot 2 \pi \cdot (k-1) \cdot (n-1)}{N}}, \quad\textrm{where}\quad 1\leq k\leq N\\
x[n] &= \frac{1}{N} \cdot \sum_{k=1}^{N} X[k] \cdot e^{\frac{+j \cdot 2 \pi \cdot (k-1) \cdot (n-1)}{N}},\quad \textrm{where}\quad 1 \leq n\leq N
\end{align}
Scaling by signal length

My peers typically scale the data by $\small \frac{1}{N}$ immediately after the processing the fft.
(We do not consider the raw fft data before scaling.)
%% Perform fft
X_f = fft(x, n_sample, 1) / n_sample;
% fft must be normalized by the number of samples in the data.
% This convention was set by the software developer (Mathworks).

Is this correct?

If so,  why does the MATLAB ifft function expect that we have not scaled by $1/N$ already?
Is there a MATLAB ifft function or function option which does not automatically scale by $1/N$?

Alternatively, is there a better convention which we should be using in placing the $1/N$? For example, placing the $1/N$ in the fft rather than the ifft, or placing an $1/\sqrt{N}$ in both equations, instead of an $1/N$?
Scaling by sampling period

I have heard that the fft and ifft functions assume that the sampling period $T_{\rm sampling} = 1/f_{\rm sampling} = 1$, and that for the functions to be true, the following would need to apply:

\begin{align}
X[k] &= \frac{1}{T_{\rm sampling}} \cdot \sum_{n=1}^{N} x[n] \cdot e^{\frac{-j \cdot 2 \pi \cdot (k-1) \cdot (n-1)}{N}},\quad \textrm{where}\quad 1 \leq k \leq N\\
x[n] &= \frac{T_{\rm sampling}}{N} \cdot \sum_{k=1}^{N} X[k] \cdot e^{\frac{+j \cdot 2 \pi \cdot (k-1) \cdot (n-1)}{N}},\quad\textrm{where}\quad 1 \leq n \leq N
\end{align}

See links:

Link 1 (see comment to Matt Szelistowski by Dr Seis)
Link 2 (see answer  by Rick Rosson vs that of Dr Seis)
Link 3 (see comment by Matt (Message: 7/16) and comment by Poorya (14/16)
Link 4 (see pg. 10, slide [1,1])
Link 5 (see pg. 8+9) [it seems he is using inverse convention for fft and ifft].

Is this true?

I'm particularly piqued because I cannot find any DFT or DTFT equations on Wikipedia which include the sampling period.

Whether or not to scale the forward FFT by 1/N depends on which result you want for further analysis: energy (preserving Parseval's identity), or amplitude (measuring height or volts, etc.).

If you want to measure or analyze energy, then don't scale by 1/N, and a longer sinusoid of the same amplitude will produce a larger FFT result, proportional to the greater energy of a longer signal.

Slightly more commonly, if you want to measure or analyze amplitudes, then to get a longer sinusoid (thus with more total energy at the exact same amplitude) to produce about the same FFT result as a shorter signal, you will need to scale down the FFT summation by an ratio proportional to the length.  The ratio could be reference_length/N, which is sometimes 1/N if the system input gain is 1.0 for whatever dimensions or units, including the time interval dimensions, that you choose to use in your further analysis.  You need to scale down proportionally because a DFT is a summation: the more you sum up similar items, the bigger the result.

So.  Energy or amplitude.  Which do you want?

Now if you scale down the forward FFT, then you should not scale the inverse so that IFFT(FFT(x)) == x.  Or do vice versa.

The 1/sqrt(N) for scaling seems to me to be for either when one needs a formal symmetry for some proof, or when building some sort of hardware pipeline where the latency and/or number of arithmetic units/gates for the DFT and for the IDFT needs to be identical.  But you get neither a good direct measurement of either energy or amplitude for any typical type of engineering analysis.


What is full AM used for? What are some applications where full AM is used nowadays? Is it used at all considering the fact that it is less power-efficient than, say, AM with suppressed carrier?

AM is still used worldwide in civil aviation, where there is a large installed base of radios (in small planes, etc.) than can not be easily updated (due to multiple national and International regulatory hurdles), and where you want one transmission to be able to talk over another signal on the same channel in an emergency.

AM radio receivers are simple and cheap using primitive technology (1 discrete diode, etc.), and there are still a few parts of the world where people can barely afford an old one of those.


Why Fourier series and transform of a square wave are different? Here is a square-wave presented by Fourier series perspective:

Above coefficients shows that a square-wave is composed of only its odd harmonics.

But here below a square-wave is presented by Fourier transform perspective:

Above plot shows that a square-wave is composed of all frequencies not only harmonics, plot is continuous.

When I look at the FFT of a square-wave it looks like the Fourier transform which is continuous.

Series and transform gives different interpretation of a square wave. Why is that?

The Fourier series expansion of a square wave is indeed the sum of sines with odd-integer multiplies of the fundamental frequency. So, responding to your comment, a 1 kHz square wave doest not include a component at 999 Hz, but only odd harmonics of 1 kHz.

The Fourier transform tells us what frequency components are present in a given signal. As the signal is periodic in this case, both the Fourier series and the Fourier transform can be calculated, and they should tell us the same information. The Fourier transform of a continuous periodic square wave is composed by impulses in every harmonic contained in the Fourier series expansion. Maybe this picture from Oppenheim's Signals and Systems may help.

The actual Fourier transform are only the impulses. The dotted-line is a sinc function that doesn't apply to this question, but gives the notion that this transform has something to do with the transform of a square pulse (i.e. a not periodic signal), which happens to be a sinc.

To put it mathematically:

The Fourier series coefficients are $$\frac{\sin(k\omega_0 T)}{k\pi}$$
The Fourier transform is $$\sum\limits_{k=-\infty}^{\infty}\frac{2\sin(k\omega_0 T)}{k}\delta(\omega - k\omega_0)$$

So the series coefficients and the Fourier transform are the same, except that there is a proportionality factor of $2\pi$ and, in the first case, you plot bars (as the coefficients do not describe a function, they are just numbers), but in the second one you have impulses (because the Fourier transform is a function).


Disjoint ROCs of input and system function I am trying to understand how the output behaves when the input and the system function do not have a common region of convergence (ROC) for an LSI system.
Consider an LSI system with $x[n]$, $h[n]$, $y[n]$ as input, impulse response, output signal respectively:
We know that
\begin{gather}
Y(z) = X(z)\cdot H(z)\\[5pt]
\text{with } \ ROC(Y) \supseteq ROC(X) \ \cap \ ROC(H)
\end{gather}
So what if $ROC(X) \ \cap \ ROC(H)=\phi$, what would be the output? Is it even properly defined?
I tried with some examples considering $x[n]$ and $h[n]$ with disjoint ROCs, the output is not defined,
here is one example:
\begin{gather}
x[n] = (1/4)^nu[-n-1] \\[5pt]
h[n] = (1/2)^nu[n]\\[5pt]
ROC(X): |z|<1/4\\[5pt]
ROC(H): |z|>1/2\\[10pt]
y[n] = \sum\limits_{k=-\infty}^{\infty} (1/4)^ku[-k-1](1/2)^{n-k}u[n-k]\\[5pt]
y[n]= (1/2)^n\sum\limits_{k=-\infty}^{-1} (1/2)^ku[n-k]\\[5pt]
y[n] = (1/2)^n\sum\limits_{k=1}^{\infty} 2^ku[n+k]\\[5pt]
\text{Clearly } y[n] \text{ is not defined}
\end{gather}
My question is whether this is true in every case that the output would not be defined or there are some special $x[n]$ and $h[n]$ with disjoint ROCs where $y[n]$ is properly defined.

If the ROCs don't overlap, the corresponding convolution sum doesn't converge, hence the output is undefined. Note that in your question it is assumed that the $\mathcal{Z}$-transforms of the input as well as of the system's impulse response both exist.
You can construct examples for which neither the $\mathcal{Z}$-transform of the input nor of the output exist, yet the input and output signals are perfectly defined. E.g., the input sequence $x[n]=A\cos(\omega_0n)$ has no $\mathcal{Z}$-transform, neither has the corresponding output sequence $y[n]=B\cos(\omega_0n+\phi)$, where $B$ and $\phi$ depend on the system's frequency response. [In this example it is assumed that the ROC of the system's transfer function $H(z)$ includes the unit circle, i.e. the system is stable.]


Logic behind unsharp filter in matlab I am trying to understand the logic behind unsharp filter in matlab. If I type:

image = imread('nameoffile');
filter = fspecial('unsharp', 0.5);
resultImage = imfilter(image, filter);

resultImage is sharper than image. I opened fspecial.m file and saw that unsharp filter is created as:

 h = [0 0 0;0 1 0;0 0 0] - fspecial('laplacian',alpha);

If I filter an image with this matrix:

[0 0 0;
 0 1 0;
 0 0 0;]

result image is the same with original one. I know laplacian filter detects sharp edges. So, if we subtract sharp edges from an image we must get blurred image, I suppose. However, if we filter an image with h matrix(unsharp filter), we get sharper image. I did not understand the logic.

The Laplacian filter does actually outputs values of 0 on edges, and positive/negative inside image regions. Edges are detected as zero-crossings of the Laplacian.
This makes sense if you remember that the Laplacian is a second-order derivative, and not a first order.

So, the unsharp filter that you describe computes (identity - Laplacian) convoluted with the image, i.e., the image minus something that is zero on edges. Hence, high frequencies (edges) are not attenuated while lower frequencies are weakened.


What exactly is a 90 degree phase shift of a digital signal in FM demodulation appraoches? I am working on an FSK demodulator (1200/2200Hz, 1200 baud) featuring a 90 degree phase shift operation. I don't exactly understand what 360 degrees mean for a digital signal.

My sampling frequency is 8 kHz. Is a 90 degree shift a delay by 8000/4 samples, or is it a delay by 4000/4 samples?

The OP asked about a phase shift specifically but from the written details without seeing the actual demodulator implementation, I suspect he may possibly be asking how to implement the delay that is used for some simple FM demodulation approaches. I want to mention the possibility that this may actually be a delay element and not a phase shift in the sense that phase shift implies a shift of that phase at all frequencies, while a true delay has a phase shift that is proportional to frequency. It is is this property that is exploited to form simple FM demodulator structures.

For example, a very simple non-coherent digital (or analog!) FM discriminator is implemented by a delay and multiply. As introduced above, this is truly a delay resulting in a phase shift of 90 degrees only at the center frequency and NOT a phase shift of 90 degrees at all frequencies within the signal bandwidth (as would be provided by the Hilbert Transform).

As the OP commented in another answer to this question, a delay of 1.176 samples (which is indeed conveniently a true delay!) would work for this approach, the details and trade space involved will be clearer once the demodulation process of converting frequency change to amplitude change is understood.

To understand the operation, first it is useful to review how a multiplication in time between two signals results in a phase measurement:

In the case of a real signal (which I assume you may be doing), the product creates the sum and difference of the two input frequencies (and phase), and once low pass filtered, the difference signal results in a voltage that is proportional to the cosine of the phase. The complex case is also shown (and in this case either the real or the imaginary output as shown could be used, or both for non-ambiguous phase from 0 to 360 degrees).

The ideal operating point for a sensitive phase to voltage (or magnitude regardless of units) for the case of a real signal is when the two input signals are in quadrature phase; at this point the output signal is crossing through zero, but importantly, has the highest rate of change to any change of phase in the input. This is particularly true for analog demodulation, where we desire a linear phase to voltage translation to minimize harmonic distortion (so not as critical in FSK demodulation).

The next important item to understand is how a fixed delay is a "frequency to phase converter". A fixed delay passes all frequencies with equal magnitude, but the phase at the output (relative to the input) is linearly proportional to the frequency (in a fixed delay a higher frequency signal may take several cycles to transition the delay while a lower frequency signal may only shift fractions of a cycle). Therefore with regards to a signal that is changing its frequency over time (the FM modulated signal including FSK), the delay and multiply will first convert frequency variation to phase variation via the delay and then phase variation to amplitude variation via the multiplier!

It can be further determined from this understanding what delay would be optimum for two FSK frequency choices that will provide a peak maximum and minimum signal at each location (refer to the figure below which shows the complex multiplier approach but applies equally to the real signal approach with the additional low pass filter as given in the figure above.)

So for the OP’s case specifically with a 1 sample delay as the delay element:

Frequency    Phase Shift for 1 sample delay

1200 Hz:      360*1.2/8 = 54°   cos(54°) = 0.59

2200 Hz:      360*2.2/8 = 99°    cos(99°) = -0.16

This will work in the sense that it provides a bipolar output ranging from -0.16 to +0.59 out of the possible range of +/-1. So not ideal in terms of best SNR but will certainly work.

Interestingly consider the case where a 4 sample delay is used for this example:

1200 Hz:      360*1.2/8 * 4 = 216°   cos(216°) = -0.8

2200 Hz:      360*2.2/8 * 4 = 396°    cos(396°) =  +0.8

This is much closer to the ideal result for this non-coherent demodulation approach of +/-1 output and very simple to implement given the sampling rate and FSK frequencies used. As an even further simplification in a mixed digital/analog approach, a single XOR gate could be used at the multiplier, followed by an RC low pass filter (and then a Schmidt trigger as your decision slicer)! It is important in all of these implementations that the input signal is hard limited (after bandpass filtering) to remove any incidental AM modulation/noise on the signal prior to demodulation, since these demodulation approaches are equally sensitive to amplitude variation on the input.

To add to this because I like the simplicity and the similar design completely in the analog domain, below shows a commonly used non-coherent FM demodulator circuit. This is not an optimized FM demodulation solution but commonly used due to it simplicity.

Observe from the sketch below the expected magnitude and phase response of just the LC resonator portion of the circuit (the capacitor and inductor in parallel and connected to ground):

At low frequencies the inductor impedance ($j\omega L$) dominates resulting in a positive 90° phase shift (+j), and at high frequencies the capacitor impedance dominates ($1/(j\omega C)$, resulting in a -90° phase shift (1/j). The important thing to observe is what the phase shift is doing as the response passes through resonance at $f_c$. Near this location the phase change versus frequency is nearly linear with a slope defined by the Q of the resonator ($\sqrt{L/C}$). The slope, or derivative of phase with regard to frequency is time delay! So the tank is creating a delay line in proximity of this location (resonance), and conveniently at resonance the impedances cancel resulting in no loss (in the ideal circuit). However unfortunately the phase shift at this frequency, as shown in the response, is 0°. This is solved by adding another smaller capacitor (smaller to dominate in impedance) in series which serves to shift the entire phase response down an additional -90°!  So the circuit gives us the desired delay line with a 90° phase shift at the center frequency (carrier) of our modulated signal. We choose the slope (Q) of the tank circuit based on the discriminator sensitivity needed, and possible linearity specs that we need to meet (as for larger frequency excursions relative to the bandwidth of the tank circuit used the phase versus frequency will no longer be linear and there will be FM to AM translation effects due to the amplitude response of the tank).


How can I find SNR, PEAQ, and ODG values by comparing two audios? I have two audios, A and B. B is the noisy version of A. I want to know how much B is different from A.
Can I do this by calculating signal to noise ratio (SNR), Objective difference grade (ODG), and Perceptual Evaluation of Audio Quality (PEAQ)?

How can I calculate them in Matlab?

I recommend calculating the correlation coefficient; the correlation coefficient will directly tell you how much B is (linearly) different from A, and you can compute the SNR in dB from the correlation coefficient $\rho$ using

$$SNR = 10log_{10}\left(\frac{\rho^2}{1-\rho^2}\right)$$

Given you have a noise-free version of your waveform, this would be a perfect application for the approach above. In this case "noise" is any deviation from your noise-free waveform once the two signals are properly scaled (which the computation will do).

See this similar post where I detail the correlation coefficient calculation:

Noise detection


Correct color values on image taken with no infrared filter I have a Raspberry Pi NoIR Camera. Which I want to take pictures with at daytime in a well lit room.

Unfortunately the colors on the pictures do not look anything like you would expect from a normal camera. I guess the reason is:

  (NoIR = No Infrared.) This means that pictures you take by daylight
  will look decidedly curious, but it gives you the ability to see in
  the dark with infrared lighting.

This post states that you can only overcome this issue during daylight by using an additional infrared filter or by using a "normal" camera.

So the pictures seem to lack any green values, plus some of the black values appear red/purple.

I could apply filters to the images programmatically or using gimp and I was wondering how would one tackle such a task? Could this be solved with some simple color correction algorithms? Is it even possible to revert this effect without using a additional "hardware"-infra red filter?

Thanks in advance

Unfortunately, no practically known signal processing method could provide what you wanted yet. Instead, you must use a hardware optical IR filter that would block the incoming IR waves, before they reach into your sensor.
An optical image is an array of intensity values captured by the sensor pixels which are illuminated by the incoming electromagnetic wave during the exposure time. During the capturing process, the static image sampling is performed on the 2D sensor surface, unlike an audio sampling which captures the time based variations of the incoming acoustic pressure wave.
Now the static-image sampling does not capture the time variations of the incoming electromagnetic wave. So those array of pixels represent only the spatial distribution of the illumination intensity $I(x,y)$ accross the sensor surface, over which an image was formed by a focusing lens system. Therefore, time-based property of light's wavelength is not represented in the samples, and cannot be processed later for modification.
Note that the camera system also samples the wave along temporal (time) dimension, the motion-picture; however it's typically about a few hundred frames per second at max. Given the required bandwidth of the visible spectrum; about $10^{14}$ Hz (visible light in [400-700] nm, plus IR spectrum to at least 1000 nm for near IR) this would require a sampling rate on the order of $10^{14}$ Hz. Such high data-rate is not practical for commercial applications of typical imaging devices.
So, for your Sony NoIR cam, the temporal-spectral information of the incoming electromagnetic wave is not captured, not represented by the image pixels, and therefore cannot be processes by any further means.
You should therefore physically place an IR blocking filter in front of the camera sensor, or the lens. Old GSM phones with built-in cameras, may incorporate an IR filter..?


Image Processing: How to detect a quadrialteral signboard in image? How can I detect a quadrilateral signboard from an image captured using a mobile phone?
How can i detect shapes like rectangle? rounded rectangle (rounded corners instead of shape corners)?

I am using opencv.wrapper  but i am new to it.

Thanks.

Here is the sample:
alt text http://www.freeimagehosting.net/uploads/b03442fd36.png

alt text http://www.freeimagehosting.net/uploads/e6b36040e8.png

Due to the noise and many lines there, i am not able to determine which is the boundary line of the signboard. Sometimes I can find the boundary of the lines after hough transform. I am trapped....in this kind of scenario...

This is the 2 raw pictures taken from mobile phone camera

alt text http://www.freeimagehosting.net/uploads/6dbd613edf.jpg
alt text http://www.freeimagehosting.net/uploads/720da20080.jpg

I need ya advice to see how i can process the image to get the signboard out?

Thank you very much

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

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

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

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

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

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

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

The result looks like this:

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

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

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

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

Results:


Finding the right measure to compare sound signals in the frequency domain I am comparing the spectral contents of some sound signals with MATLAB. I have quite a few measures but wanted one that represents the centre of gravity of the spectrum.

I used SpecCentroid() found in File Exchange to calculate the linear mean spectrum centroid. Then I also used the function meanfreq() and medianfreq() from MATLAB.

I am struggling to understand (mathematically and perceptually) the differences of these measures and so I am not too sure which one to choose to best compare signals between each other.

Thanks a lot.

EDIT:

I have 2 types of sounds from different sources and try to compare their spectral 'signatures' to predict if they both would be relevant to some animal species. At this point, I am only interested in the frequency domain (although ultimately the time domain will be as equally important). I have so far compared the peak frequencies. However, as the signals are quite noisy, I was wondering if other measures could be more appropriate. That is why I was exploring the centroid, mean, median options.

Here is the same plot as before in the linear scale:

EDIT 2:

After Laurent Duval's answer and suggestions, I tested the estimated peak location, with different weights. For the signal example given above, I got:

\begin{array}{|c|c|c|}
\hline\text{$p$} & \text{$f$} & \text{dispersion $d$} \\
\hline
1 & 979.6 & 979 \\
2 & 493.12 & 627\\
3 & 327 & 295\\
4 & 296 & 135\\
5 & 291.14 & 65\\
6 & 285 & 35\\\hline
\end{array}

The peak frequency as found with findpeaks() was $285$.

Let's take one of my more noisy signals:

\begin{array}{|c|c|c|}
\hline\text{$p$} & \text{$f$} & \text{dispersion $d$} \\
\hline
1 & 1460 & 826\\
2 & 1295 & 768\\
3 & 1122 & 719\\
4 & 955 & 648\\
5 & 811 & 562\\
6 & 698 & 470\\\hline
\end{array}

The peak frequency as found with findpeaks() was $427$.
Example at $p=4$:

It is quite clear that the dispersion is much higher than for the last (less noisy) signal.

[EDIT: a related discussion on (power) spectrum centroids can be found in How to calculate the mean/center frequency of the spectrum?]

Let us assume you have a relatively concentrated peak, with noise. I assume your maximum frequency is ok. If too high, you can trim the spectrum to lower frequency, which could help.

The code below gives a simple model, estimates a peak location, and a gross measure of dispersion (or uncertainty). It is based on a weighted average or center of mass. If the peak is narrow, the center of mass (and the median) will be drifted away toward the center of the frequency axis ($f_i$'s). Because even if low, the noise spreads. This is an illustration of Archimedes' lever principle:

  Give me a place to stand on, and I can move the earth.

To achieve your goal, you can give more weight to your spectrum of amplitude $a_i$, and compute an estimated peak location $\bar{f}$:
$$\bar{f} = \frac{\sum_i w_i f_i}{\sum_i w_i}\,. $$

If you choose $w_i = a_i$, you get the standard average. This is not enough. You can increase the peak weight with a power: $w_i = a_i^p$. This is the aim of weightPower = 4 in the code. The location is given by the red o, and the x are the left and right bounds derived from a weighted standard deviation. The higher the $p$, the sharper the estimation, with a risk of major errors if the data differs from the model.

Please report your tests, to adjust the answer. And the result follows. With $p=1$, you get $360$, with $p=2$, you get $243$, with  $p=4$, $\bar{f}=201$, quite close to $200$ in my model. If you need much higher robustness, you can go for weighed medians instead of weighed means.

freqAxis = linspace(0,1000,1000)';
freqPeak = 200;
freqPeakAmplitude = 1;
freqGGDExponent = 2;
freqGGDWidth = 50;
freqNoise = 0.2;

dataSpectrum = exp(-((freqAxis-freqPeak)/freqGGDWidth).^freqGGDExponent)+freqNoise*rand(size(freqAxis));

weightPower = 4;
freqEstLocation = sum((dataSpectrum.^weightPower).*freqAxis)/sum(dataSpectrum.^weightPower);
freqEstDispersion = sqrt(sum((dataSpectrum.^weightPower).*(freqAxis-freqEstLocation).^2)/sum(dataSpectrum.^weightPower));

freqEstAmplitude = interp1(freqAxis,dataSpectrum,freqEstLocation);
clf;hold on
plot(freqAxis,dataSpectrum);
plot(freqEstLocation,freqEstAmplitude,'or');
plot([freqEstLocation-freqEstDispersion freqEstLocation+freqEstDispersion],freqEstAmplitude,'xr');
axis tight;grid on


Is it possible to define Fourier or wavelet transforms on DNA sequences? I am wondering how and if it is possible to define a Fourier transform or Wavelet transform on DNA sequences which are basically arrays with the values $\{T,C,G,A\}$ in them.
I have found a paper which uses a "Voss 4D binary indicator representation":

Yin, C., Chen, Y. and Yau, S.S.T., 2014. A measure of DNA sequence
similarity by Fourier Transform with applications on hierarchical
clustering. Journal of theoretical biology, 359, pp.18-28.

Their code is supposed to be available online FFTDNA4D.m, but the link is broken.
The general idea is to define an indicator mapping per nucleotide type ($\alpha \in \{T,C,G,A\}$) for the DNA sequence $s(0), s(1), \dots, s(N-1)$
$$
\begin{equation}
  u_{\alpha}(n)=\begin{cases}
    1, & \text{if $s(n)=\alpha$}\\
    0, & \text{otherwise}
  \end{cases}
\end{equation}
$$
Then the DFT is given by:
$$
U_x(k) = \sum_{n=0}^{N-1} u_x(n) e^{-2\pi i \frac{kn}{N}}
$$
the power spectrum is then given by:
$$
\Phi(k) = \sum_{x\in \{ A, T, G, C\}} |U_x(k)|^2 \quad k = 0,1, \dots, N-1
$$
Is there a better way to do it? This method masks out the differences between the different types of DNA bases.

First, a code for paper A measure of DNA sequence similarity by Fourier Transform with applications on hierarchical clustering. Journal of theoretical biology is available at GenomeDFT (github).
Then, on the transforms. Discrete Fourier transforms can be defined on other domains than vector spaces. There are Fourier transforms on groups, rings or finite fields. On the latter, it is commonly called a number-theoretic transform (NTT). This possibility entails to making sense of group/ring/field axioms for the nucleotides. For instance, it should be possible to use a finite field of characteristic four. Apparently weird mathematical structures could be interesting when considering (C,G) or (A,T)  patterns, or the redundancy in the codon/amino-acid correspondence.
Some papers mention the Ramanujan-Fourier transform (A Novel Method for Comparative Analysis of DNA Sequences by Ramanujan-Fourier Transform, apparently with the same first author of the paper you mentioned). With less structured transformations, it is also possible to evaluate periodicity, as for instance in Categorical spectral analysis of periodicity in nucleosomal DNA.
There also were some works on wavelets, either continuous or discrete. On the first flavor, by A. Arneodo and colleagues, and more recently in 2019 with
Wavelet-Based Genomic Signal Processing for Centromere Identification and Hypothesis Generation. With discrete wavelet transform, there is for instance Identification of exonic regions in DNA sequences using cross-correlation and noise suppression by discrete wavelet transform, Application of discrete wavelet transform for analysis of genomic sequences of Mycobacterium tuberculosis, Wavelet analysis of DNA sequences.
Finally, on the mappings. There is for instance Digital signal processing methods for biosequence comparison, 1991 or Mapping Equivalence for Symbolic Sequences: Theory and Applications, 2009.
As asked by Marcus, your actual purpose might direct you to an appropriate combination.


Gaussian filter as a low pass filter Can anyone explain me how is a Gaussian filter a low pass filter? It may be a simple thing but I just can't seem to wrap my head around it. Also while applying a low pass filter for bandlimiting (to prevent aliasing), which of the following two is better:
(a) Applying Gaussian filter to the signal.
(b) Applying Box filter in the frequency domain and then
      using IFT on the resultant.
I think option (a) is better because we just convolve the signal with a Gaussian while in option (b) we first need to get the FT, then apply the filter and then convert back using IFT. But what about the quality of the resultant signal?

I'm still a newbie so answers with less jargon will be appreciated.

it's because all of the filter coefficients (which is the impulse response $h[n]$) are positive and there is no "ringing" in the impulse response.

so, assuming the impulse response is scaled so that all of the $h[n]$ add to 1, the output of the gaussian filter will be a weighted average of the input.  averaging is a process that removes high-frequency variation.  so what is left are the low frequencies, like DC.  that's what a low-pass filter does.


Periodicity of complex exponential in continuous and discrete time (Eq 1.51, Signals and Systems by Oppenheim & Wilsky) Hi All: This is very basic but I've always wondered about it and now I see it in print in a textbook so I may as well ask. In Signals and Systems on page 26, it says
$$e^{j(\omega_0 + 2\pi)n} = e^{j2\pi n} e^{j\omega_0 n} = e^{j\omega_0 n}     \tag{1.51} $$
Then it says:

From eq. (1.51), we see that that the exponential at frequency $\omega_{0} + 2\pi$ is the same as that at frequency $\omega_{0}$. Thus, we have a very different situation from the continuous time case in which the signals $e^{j\omega_{0} t}$ are all distinct for distinct values of $\omega_{0}$. In discrete time, these signals are not distinct, as the signal with frequency $\omega_{0}$ is identical to the frequencies with $\omega_{0} \pm 2\pi$, $\omega_{0} \pm 4\pi$
and so on.

The reason why the concept above has always confused me is because if we replace $n$ with $t$, then it is true for integer values of $t$ also. Why does this not matter ? I know  that  concept also has to do with why aliasing cannot exist in the continuous case so, if I can understand this, a lot of things would get clear for me. Thanks for your help and the great answers in general from many.

The simple answer is that
$$e^{j\omega_0t}\neq e^{j(\omega_0+2\pi)t}\tag{1},\qquad t\notin\mathbb{Z}$$
Since $t$ is a real variable, the inequality is true for uncountably many values of $t$. Equality is only achieved for countably many integer values of $t$. From $(1)$ it follows that for real $t$, the function $e^{j\omega_0t}$ is not $2\pi$-periodic in $\omega_0$, unlike in the discrete-time case.
Two functions of a real variable can be equal on a countably infinite set of values of the independent variable, even though they're not the same function.


Computer Vision - Gaussian Pyramid, Scale Space and Pattern Matching I am planning to implement a pattern matching algorithm using something like correlation as a matching metric.

I know that the template I am going to use will, if present, have different sizes in the target images. The interval of sizes might be something like :
x 0.5 , x 0.75 , x 1.0 , x 1.25 , x 1.5 , x 2.0 of the original template size.

To speed up the search it is suggested in various articles to build a gaussian pyramid of both the template and the target search image.

However I am not able to relate in any meaningful way the number of pyramid levels to create (for model and for target) , the object scales as described above, and the gaussian filter sigmas (used when creating the pyramid) !

Could someone shed any light on this ?

For example 1 level of the pyramid reduces the resolution of the image by 2 (is this the same as saying that the image has been scaled by 0.5 ) ?

Than what about a pyramid level of -1 (which should make the image bigger x 2.0) : with what gaussian filter should it be treated before upsampling ?

What about object scales (e.g. x 1.25) which fall in-between pyramid levels ?

Thanks to anyone who could provide some insight/references,
    Todor

Template matching is related to finding the maximum correlation between two signals, where one of the signals (the image being explored) is usually much bigger than the reference signal (the object being searched for). Obviously, the best correlation score can only be obtained when the size of the object in the image is equal to its size in the template (intuitively, both shapes will superpose).

  However I am not able to relate in any meaningful way the number of pyramid levels to create (for model and for target)

The number of scales is related to two things:

the size of the object of interest in the explored image
the computation time reduction that you want to obtain.

To have the best localization accuracy, you need to have one template image that matches exactly the size of the object in the explored image. Thus, if you don't know a priori this size, but you have instead a range of sizes, your template pyramid should include each of these sizes.

Since the template matching part can be expensive (one correlation computation for each possible position of the template in the image), you want to reduce the number of pixels explored. Thus, you create a pyramid of the image to explore, start by looking for the object at the coarsest scale to get a rough position that gets refined at the next scale, etc. Again, since the size of the object needs to match between the template and the explored image, you need to rescale appropriately each of the images from the pyramid that you built for the template (previous paragraph).

  For example 1 level of the pyramid reduces the resolution of the image by 2 (is this the same as saying that the image has been scaled by 0.5 ) ?

It depends... Reducing the resolution by two along each dimension divides the number of pixels by 4. This is usually referred to as scaling by 0.5 but as an image resolution reduction by 4.

A common practice is to divide the number of pixels by 2, yielding a $\sqrt{2}$ factor along each dimension.

  Than what about a pyramid level of -1 (which should make the image bigger x 2.0) : with what gaussian filter should it be treated before upsampling ?

You need to use image interpolation techniques here. OpenCV uses (in optical flow methods) a gaussian upsampling by injecting $0$'s along the odd lines/columns then applying a gaussian blur.


Why is $\delta[an]=\frac 1 a \delta[n]$ in discrete time? Why is $\delta[an]=\frac 1 a \delta[n]$ in discrete time? Prove.

Hi, I want to prove it but I don’t know on what facts to rely. In continuous time we have integral and dt in that integral that turns to be dt/a for continuous time delta function. But in discrete time, we don’t have dt. So how to prove it then?

It's actually good that you couldn't prove it because the claim is wrong. Note that in discrete time we have $\delta[n]=1$ for $n=0$ and $\delta[n]=0$ for $n\neq 0$. So if $a\neq 0$, you simply get $\delta[an]=\delta[n]$ (note that $a$ must be an integer for this to make sense). If $a=0$ then $\delta[an]=1$ (for all $n$).


3D wiggle plot for an analytic signal: Heyser corkscrew/spiral Just reading The Analytic Impulse, A. Duncan, 1988, I met the name "Heyser corkscrew" for the first time in my DSP life, for a 3D display of a cisoid or complex exponential $e^{i\omega }$ (often denoted  $\text{cis}\, x$ for "cosine plus $i$ sine").

The view is mentioned as Euler formula's 3D view in the context of circular polarization:

This representation is quite common for analytic signals, but this name does not seem standard. Questions are:

What is the history of this Heyser name, with a precise reference?
What are the (earlier) history and names for this representation  (going back to the Fourier oblique view)?

After @MBaz hints, the term "Heyser spiral" seems a little more common. It is named after Richard C. Heyser, and can be interpreted as a 3D variant of Nyquist plots... Indeed, in  A Matter Of Frequency: The Nyquist Plot Explained, one can read:

  Dennis Gabor proposed the “analytic signal” which Richard Heyser
  developed into what has since been named the Heyser Spiral. The
  Nyquist plot is the “end view” shadow of the complex analytic signal

Looking at ancient geophysical papers, I have found in Truncations and phase relationships of sinusoids, Philip L. Jackson, Journal of Geophysical Research, 1967, this "Fourier transform" with "oblique view", that seems to predate Heyser corkscrew/spiral.


Match an Equalizer Curve I’ve passed a sine sweep wave file into an equalizer plugin I have, now I want to model it in matlab so that I can ‘copy’ the equalizer setting sound. I’ve both in.wav and out.wav audio files and Matlab but I don’t know how to extract a transfer function and then build a model to exactly recreate it using iir filters. Can someone help me please?

The OP is trying to do what we would refer to as channel estimation. I am linking two posts below where this is all worked out in greater detail but do want to make the point first that the choice of "sounding pattern" is important. The channel can only be estimated where there is signal energy in the frequency domain, so the sine sweep is an excellent choice as long as the sweep is over all the frequencies of interest. Other sounding patterns commonly used are pseudo-random noise generated with PN sequences (from linear feedback shift registers: LFSR).
An optimized FFT chirp for channel estimation:
How can I plot the frequency response on a bode diagram with Fast Fourier Transform?
Using Wiener-Hopf equations for least-squares channel estimation:
Compensating Loudspeaker frequency response in an audio signal
Related posts on PN sequences as an alternate sounding pattern:
What are the advantages and disadvantages of PN sequence over Walsh code
Sequence Length for a Linear Feedback Shift Register


Detecting defects/differences in printed images I've been experimenting with various methods to detect defects in printed images. Features that are qualified as defects may include scratches, discoloration, etc.

I'm not very experienced in image processing and not sure if my approach would be acceptable.

There are several approaches that I've tried:

Given two images (one good, the other one is defect)

Threshold and blur
Run Canny edge detection
Remove noises using fastNlMeansDenoising
XOR the two images and use SIFT to find outstanding features that may be classified as a defect.

It works very well - in very limited cases.

Seems like this method is used/proposed in industry a bit (?) as seen in this paper. But the problem with it is that if you have two images that are misaligned to the slightest bit, then you would get a problem because all your XOR result in an XOR-ed image that is very wrong. So to me it doesn't seem like an acceptable solution, since it is impossible to perfectly align two images every time (or I could try to realign the images).

How should I approach this. Is there a method that is similar to XOR-ing that is more tolerant to subtle differences? Or perhaps, a totally different approach.

Possible alternative

Break up both images into chunks of n by n pixels. Compute and compare the histogram of each chunk of pixels broken up from image A to image B. Even if one of the images is blurred a little, the histogram should be still be similar (not sure if this is a sound assumption). If any chunk contains dark or white spots, this should flag that chunk for further inspection.

Detects darker font color

Detects eyebrows that aren't supposed to be there

Detects faded/lighten area

All the comments in the comment section were very constructive, thank you.

In the end, I ended up using the solution here

For my case, most of my images are not expected to tilt more than a few degree, lighting and depths are strictly controlled. So a Euclidean Motion model (MOTION_EUCLIDEAN) should suffice for findTransformECC which computed the transformation matrix that is required to correct the misalignment.


BER based transmitter gain Are there any direct relation formula between BER and Transmitter gain. The higher the transmitter gain is, the better the BER performance is. In here, I want to plot graph for the BER based on different transmitter gains and compare the MATLAB simulation vs experimental results.

To simplify, let's consider the transmission of a single pulse $$s(t) = A p(t),$$ where $p(t)$ has energy 1, and the transmitted energy is $A^2$.

Note that the value of $A$ is related to the transmitted energy. If the transmitter's amplifiers have a lot of gain, then $A$ will be large. The AWGN channel model does not make explicit mention of this, but you should always keep it in the back of your mind that it's only a model, in a real system there is a lot more going on.

Now, the receiver has a filter matched to $p(t)$; when sampled at the right time, the matched filter output is equal to $A$. The receiver's noise has PSD equal to $N_0/2$, which is constant and independent of $A$.

You can then easily conclude that increasing the transmitted power will increase the SNR at the matched filter's output, which will in turn decrease the probability of error.

So, one possible approach would be:

Define a fixed transmitted constellation. Let's say you use BPSK with constellation $\lbrace -1, 1 \rbrace$.
Define a transmitter gain G.
Calculate the SNR at the filter's output: $2G^2/N_0$.
Plot the BER as a function of this SNR.


Advise on detecting anatomical landmarks in a CT reconstructed volume I'm trying to automatically detect some medical defined anatomical landmarks in a CT reconstructed volume. Medical doctors use these landmarks to measure some patient specific parameters.
I have attempted to use the SIFT feature descriptor, since these anatomical landmarks are kind of "keypoints". This did not work very well since the landmarks are points (or tiny regions) that are in general not "interest points" as defined by SIFT.
I have been looking many pattern/template matching algorithms but, when I do not have rotation/translation/scale problems, I find that the extracted features do not differentiate each landmark enough (from the rest of the landmarks and from the rest of the non landmark patches) to train a classifier that performs well enough (at least an 80% of detection accuracy).

Please let me know if I'm not stating the problem clearly enough.

I would really appreciate any advise.

Thanks!

Example image:

The small x crosses and little squares are over the landmarks I want to detect (I forgot to mention that I have a training set, with the labeled landmarks). The white lines represent the measures taken. These are some slices of different cases (of course, I cannot post the full 3D volume).

I hesitate to write this as an answer, but given that you're asking only for advice, I will do so.

I suggest investigating techniques based on the Dual-Tree Complex Wavelet Transform (DTCWT). These have shown to be useful for generating descriptors that have good tolerance to shift, scale and rotation of the source images. It's not the classic problem in that you're not allowing the points to be assigned for you, but I suspect with some thought you can adapt the techniques to predefined landmarks.

Clearly, the landmarks have some interest from the perspective of a clinician, so there is something of interest about them - it's simply a case of modelling that in the descriptor. Wavelet techniques (in particular the DTCWT) tend to be good at modelling features that the eye picks up on.

A starting point would probably be this fairly recent paper.


Confusions regarding convolution? As far as i have studied and understood, convolution is the process by which we can get/determine output of LTI systems
While reading one web link about convolution, i came across certain notation ,that i couldn't understand as shown highlighted in attached photo

1)Are these notations referring to delay/shift in input and impulse response?

2)Also it mentions that **Periodic or circular convolution is also called as fast convolution as shown highlighted in last line of 2nd photo attached.  Is it idea correct?**

1)Are these notations referring to delay/shift in input and impulse response?

Yes

  2)Also it mentions that **Periodic or circular convolution is also called as fast convolution

That's a bit of a misrepresentation. For a non-trivial length of signal and/or impulse response the fastest way to implement convolution is to multiply in the frequency domain. That does indeed implement circular convolution. That's typically NOT what you actually want. If you want to leverage frequency domain multiplication for linear convolution, you need to add the right amount of zero padding and/or overlap handling to your algorithm. Read up on "Overlap Add" or "Overlap Save" algorithm.


Trouble understanding an article about MFCC I'm willing to apply the same audio preparation that this article page 3.

However, I don't understand the final part :

The final 28-D output feature vector consists of the 13-D vector output from the DCT plus
the log mean energy to account for volume, and their first temporal
derivatives

Can someone tell me what's exactly in this final vector and why did the authors chose to use such a method ?

So the output of the Discrete Cosine Transform gives a 13-vector. To this, they append the log mean energy. That gives a 14-vector.

They then subtract from this the previous 14-vector to yield another 14-vector, which they concatenate with the existing 14-vector.

That gives a 28-vector.

I suspect the log mean energy is just:
$$
{\tt LME} = \log\left(\frac{1}{N} \sum_{N} x_n^2 \right)
$$
where $N$ is the length of the frame used.


Can you explain likelihood function in the equation? In the paper here, just before equation 12, authors introduce likelihood function $\pi(r|\theta)$. What does $\theta$ represent here ?. What is the intuitive explanation behind equation 12?

To me, it looks like $\theta$ is just a throw-away "any variable". They might just as well have written $r(s|\cdot)$.

Equation 12 is just saying that because $r$ and $s$ are mutually independent, the distribution of the error $\rm\bf e$ does not change. Compare (12) with (10) rearranged to give an expression for $\rm\bf e$.


IQ sampling vs 2X f0 sampling for SDR Are there any modulation schemes (or other situations) where sampling at or very near twice the carrier frequency would be better than IQ or quadrature sampling using a local oscillator at or near the carrier frequency for SDR demodulation?

I can envision advantages (or convenience) in doing this for lower IF signals certainly, as long as the sampling rate is significantly higher than the signal bandwidth to contain the phase dispersion you would get with this approach.  An implementation can be lower cost and simpler in the analog with a single mixer/ADC front-end for example. This becomes less feasible at higher carriers where nonlinear noise contributions can start to dominate (relative SFDR in ADC's for example).

To further comment on cautions with an approach of sampling at "twice the carrier frequency". I assume in this statement that the approach is to approximate the baseband IQ samples by demultiplexing the samples: I1 Q1 I2 Q2 I3 Q3, etc.  It should be pointed out that this is just an approximation, with large errors as the signal bandwidth becomes large relative to the entire digital bandwidth. This is because in this approach the 90° shift is being emulated by a delay, but a delay would only be 90° at ONLY the center frequency of the signal (and varies linearly with frequency everywhere else). To do IQ sampling without delay distortion requires a Hilbert transform which is 90 degrees at all frequencies within the bandwidth of the signal. Note this is of minimal consequence when the bandwidth is significantly smaller than the sampling rate, which would apply to most direct carrier sampling schemes, but is mentioned since the real application to do 2x carrier sampling would not be such cases but lower IF carriers where analog cost and performance is less affected.

So that said, I would likely never sample at twice the carrier frequency, but if the intention is to sample close to that (for some reason), I would choose either 2/.75 = 2.66x or 2/1.25 = 1.6x. Let me explain why below:

So how to do something close to the intention of 2x carrier sampling (which is to sample at a different frequency closest to fs with a single ADC and avoid having to do any IQ in the analog). Obviously we can envision approaches of sampling with a much higher rate to a digital IF and processing IQ in the digital, but to quantify a specific approach that has no additional inaccuracies beyond the analog hardware is to to the following:

Sample the analog signal following all rules of Nyquist (including anti-alias bandpass filtering), meaning ensure the sampling rate is twice the bandwidth of the anti-alias filter, with the additional conditions:

Sample the analog signal such that the resulting digital signal resides at $F_s/4$ where $F_s$ is the sampling rate or any integer fraction of the actual sampling rate.

This means the signal can be centered in the analog prior to sampling at $nF_s/4$ where n = 1,3,5,...    Which n is used is referred to as which "Nyquist Zone". n=1 is the first and primary Nyquist zone, for example.

Once digitized, the signal can be complex down-converted to baseband if complex IQ samples are still desired.  What is great about positioning the signal at $F_s/4$ is that the complex downconverter reduces to multiplying the signal by the following to generate I and Q:

Since the complex rotator to do $Fs/4$ downconversion to baseband is multiplying by $e^{-jn\pi/4}$ = +1, -j, -1,+j,... then the I samples for this multiplier are +1,0,-1,0,... and the Q samples are 0, -1, 0, +1 ... So baseband I and Q is created by doing the following to the sampled sequence:

$$I = +1x_1 + 0x_2 -1x_3 +0x_4 + 1x_5 ...$$

$$Q = 0x_1 - 1x_2 +0x_3 +1x_4 + 0x_5 ...$$

Very convenient: demultiplex and multiply each by +/-1!

Note, if even Nyquist zones are chosen for the analog downconversion, spectral inversion will result. This is easily resolved by rotating the other way: $e^{+jn\pi/4}$.

$$I = 1x_1 + 0x_2 -1x_3 +0x_4 + 1x_5 ...$$

$$Q = 0x_1 + 1x_2 -0x_3 -1x_4 + 0x_5 ...$$

Note how this is done just by changing the sign of the Q.

Now if the complex rotator is used for carrier tracking (which is entirely likely) then we would not use the reduction down to +/-1's in the multiplier, but still benefit by the ideal placement of fs/4 which optimizes the analog anti-alias filter design by centering the signal within the bandwidth from the edges where aliasing would otherwise occur.


Low Pass Filter in AM Demodulation I am trying to recover my baseband signal that was AM modulated. I am using the filter command in MATLAB to perform the low pass filtering ...the code is given below:

Fsig=10;                    Fc=100;              Fsamp=10.*Fc;
t = 0:1/Fsamp:2*pi;

Tx = cos(2.*pi.*Fsig.*t);               % Transmitted Signal
Carrier = cos(600.*t);                  % Carrier Signal
Modulated = Tx.*Carrier;                % Modulated Signal
Demodulated = Modulated.*Carrier;       % Demodulated Signal

Rx= filter([1 1],[1 -0.8],Demodulated);

figure(6);
plot(t,Tx,'r',t,Modulated,'b',t,Demodulated,'g',t,Rx,'k','LineWidth',1.5); grid on;
legend('Trans.','Modulated','Demodulated','Baseband')

How can I select the coefficients of the filter to perfectly recover my baseband signal?
Is there any other technique/command that can be used in matlab to recover the baseband signal other than this?

MATLAB has amdemod (see MATLAB documentation) which can be used to recover suppressed carrier AM modulated signal. From the documentation, the demodulator uses a low-pass filter generated using [num,den] = butter(5,Fc*2/Fs). This can be a good start.


How do I compute the DC component given the other part of the Fourier spectrum? I often see examples in textbooks where someone constructs a 2D Fourier spectrum of an image by fixing two (symmetric) pixels to a given value, and computes the inverse DFT to get a nice waveform image. I want to know how these images are created. It seems to me that they pick one pixel that would give them the desired spacing and orientation, fill in the symmetric opposite pixel with the conjugate value, and then somehow fill in the DC component.

This last part is what's confusing me. I want to answer the more general question: if you're given the non-DC frequency components of the Fourier spectrum of an image, how do you compute the DC component? I would be happy with an algorithmic solution, as I'm working on signals in Python with the scipy/numpy suite.

Not at all. By definition, the components are orthogonal.


How to determine if a system is LTI? If the difference equation that defines the DT system is,

$$z[n] = (1-x[n]) + (1-y[n]) $$

How can we determine in the system is LTI or not when:

$x[n]$ is the input
$y[n]$ is the input
both are inputs

Using the definitions is somewhat vague in this case, so I am lost how to explicitly explain that they are not linear and time invariant for 1 and 2 at least.

HINT:

If $z_1[n]$ is the response to input $x_1[n]$ (and $y[n]$ is NOT an input signal), and $z_2[n]$ is the response to $x_2[n]$, then it is straightforward to show that the response to a linear combination $ax_1[n]+bx_2[n]$ is NOT equal to $az_1[n]+bz_2[n]$, which would be a requirement for linearity. And due to symmetry, the same is the case if $y[n]$ is the input (and $x[n]$ is NOT an input signal). The test for time-invariance also fails in both of these cases, because only the input signal is shifted (and not the other nuisance signal).

When both $x[n]$ and $y[n]$ are input signals, the linearity test above still fails, so the system is non-linear (because it is an affine mapping due to the constants in the system equation). It is, however, time-invariant because if both input signals are shifted by some constant, the output is shifted by the same constant (which is not the case if only one of the two signals is an input signal).


What does it mean for a function to have frequencies? In a lecture, my professor mentioned that $\cos$ has two frequencies. I see that using the inverse Euler's formula we can express $\cos$ as a some linear combination of complex exponentials, each with a $\omega$ term.  But, I don't really have an intuitive understanding of what it means for a function to "have" a frequency. What other functions have a frequency, and why?

I understand that "to have frequencies" can be misleading. This is a shorthand for "having non-zero energy or amplitude at two (specific) frequencies". Rephrased, for the cosine: when a  cosine signal is represented in the frequency domain, it has two non-zero components at frequencies $\pm \omega$, and the other values anywhere on the frequency axis have zero amplitude. Or: its frequency spectrum has only two non-zero values.

Frequency in signal processing often refers to a domain of analysis well-suited to certain signals and systems (linear systems for instance). With Fourier analysis, such signals can be represented as a weighted sum of sines or cosines, and essentially the frequencies a signal has is directly  related to the number of sines/cosines with non-zero weight in the sum. One could be more precise, depending on which type of signal and Fourier transformation one is talking about (discrete/continuous).

The concept of frequency is somehow independent on the signal, just like the time, and time and frequency are interlinked via Fourier transforms. So "having frequencies" by itself has no meaning for a signal, just as "having time". But a signal can be described with either amplitudes in the time domain, or amplitudes in a frequency domain.

In that sense, whenever a Fourier transform is properly defined for a signal (integrability, etc.), and there is at least one non-zero value, the signal "has frequencies".

So, most signals you will meet have frequencies, except the zero signal, and those whose Fourier transform cannot be defined correctly. As you might see later in signal processing lectures, one often uses "having frequencies" in a loose sense, meaning "above a certain threshold" instead of being simply non-zero, as in DSP practice, it is rare to have exactly zero values (finite signals, approximations, etc.)


What Is the Relationship Between a Kalman Filter and Polynomial Regression? What is the relationship, if any, between Kalman filtering and (repeated, if necessary) least squares polynomial regression?

1. There is a Difference in terms of optimality criteria
Kalman filter is a Linear estimator. It is a linear optimal estimator - i.e. infers model parameters of interest from indirect, inaccurate and uncertain observations.
But optimal in what sense? If all noise is Gaussian, the Kalman filter minimizes the mean square error of the estimated parameters. This means, that when underlying noise is NOT Gaussian the promise no longer holds. In case of nonlinear dynamics, it is well-known that the problem of state estimation becomes difficult. In this context, no filtering scheme clearly outperforms all other strategies. In such case, Non-linear estimators may be better if they can better model the system with additional information. [See Ref 1-2]
Polynomial regression is a form of linear regression in which the relationship between the independent variable x and the dependent variable y is modeled as an nth order polynomial.
$$ Y = a_0 + a_1x + a_2x^2 + \epsilon $$
Note that, while polynomial regression fits a nonlinear model to the data, these models are all linear from the point of view of estimation, since the regression function is linear in terms of the unknown parameters $a_0, a_1, a_2$. If we treat $x, x^2$ as different variables, polynomial regression can also be treated as multiple linear regression.
Polynomial regression models are usually fit using the method of least squares. In the least squares method also, we minimize the mean squared error. The least-squares method minimizes the variance of the unbiased estimators of the coefficients, under the conditions of the Gauss–Markov theorem. This theorem, states that ordinary least squares (OLS) or linear least squares is the Best Linear Unbaised Estimator (BLUE) under following conditions:
a. when errors have expectation  zero i.e. $E(e_i) = 0 $
b. have equal variances i.e.  $ Variance(e_i) = \sigma^2 < \infty $
c. and errors are uncorrelated i.e. $ cov(e_i,e_j) = 0 $

NOTE: that here, errors don't have to be Gaussian nor need to be
IID. It only needs to be uncorrelated.

2. Kalman Filter is an evolution of estimators from least square
In 1970, H. W. Sorenson published an IEEE Spectrum article titled "Least-squares estimation: from Gauss to Kalman." [See Ref 3.] This is a seminal paper that provides great insight about how Gauss' original idea of least squares to today's modern estimators like Kalman.
Gauss' work not only introduced the least square framework but it was actually one of the earliest work that used a probabilistic view. While least squares evolved in the form of various regression methods, there was another critical work that brought filter theory to be used as an estimator.
The theory of ﬁltering to be used for stationary time series estimation was constructed by Norbert Wiener during 1940s (during WW-II) and published in 1949 which is now known as Wiener filter. The work was done much earlier, but was classiﬁed until well after World War II). The discrete-time equivalent of Wiener's work was derived independently by Kolmogorov and published in 1941. Hence the theory is often called the Wiener-Kolmogorov filtering theory.
Traditionally filters are designed for the desired frequency response. However, in case of Wiener filter, it reduces the amount of noise present in a signal by comparison with an estimation of the desired noiseless signal. Weiner filter is actually an estimator. In an important paper, however, Levinson (1947) [See Ref 6] showed that in discrete time, the entire theory could be reduced to least squares and so was mathematically very simple. See Ref 4
Thus, we can see that Weiner's work gave a new approach for estimation problem;  an evolution from using least squares to another well-established filter theory.
However, the critical limitation is that Wiener filter assumes the inputs are stationary.  We can say that Kalman filter is a next step in the evolution which drops the stationary criteria. In Kalman filter, state space model can dynamically be adapted to deal with non-stationary nature of signal or system.
The Kalman filters are based on linear dynamic systems in discrete time domain. Hence it is capable of dealing with potentially time varying signal as opposed to Wiener. As the Sorenson's paper draws parallel between Gauss' least squares and Kalman filter as

...therefore, one sees that the basic assumption of Gauss and Kalman
are identical except that later allows the state to change from one
time to next. The difference introduces a non-trivial modification to
Gauss' problem but one that can be treated within the least squares
framework.

3. They are same as far as causality direction of prediction is concerned; besides implementation efficiency
Sometimes it is perceived that Kalman filter is used for prediction of future events based on past data where as regression or least squares does smoothing within end to end points. This is not really true. Readers should note that both the estimators (and almost all estimators you can think of) can do either job. You can apply Kalman filter to apply Kalman smoothing.
Similarly, regression based models can also be used for prediction. Given the training vector, $X_t$ and you applied $Y_t$ and discovered the model parameters $α_0 ... a_K$ now for another sample $X_k$ we can extrapolate $Y_K$ based on the model.
Hence, both methods can be used in the form of smoothing or fitting (non-causal) as well as for future predictions (causal case). However, the critical difference is the implementation which is significant. In case of polynomial regression - with entire process needs to get repeated and hence, while it may be possible to implement causal estimation but it might be computationally expensive. [While, I am sure there must be some research by now to make things iterative].
On the other hand, Kalman filter is inherently recursive. Hence, using it for prediction for future only using on past data will be very efficient.
Here is another good presentation that compares several methods: Ref 5
References

Best Introduction to Kalman Filter - Dan Simon Kalman Filtering Embedded Systems Programming JUNE 2001 page 72

Presentation: Lindsay Kleeman Understanding and Applying Kalman Filtering

H. W. Sorenson Least-squares estimation: from Gauss to Kalman IEEE Spectrum, July 1970. pp 63-68.

Lecture Note MIT Course ware - Inference from Data and Models (12.864) - Wiener and Kalman Filters

Presentation Simo Särkkä From Linear Regression to Kalman
Filter and Beyond Helsinki University of Technology

Levinson, N. (1947). "The Wiener RMS error criterion in filter design and prediction." J. Math. Phys., v. 25, pp. 261–278.


Extract approximation and detail coefficients I am a student and would like to know how to extract the approximation and detail coefficients (A5, D5, D4, D3, D2 and D1) after an applied DWT decomposition dB5 (Daubechies 5, level = 5)?

This sounds quite doable, reading the online documentation, for instance:
Multilevel DWT, IDWT and SWT, by adapting the example given:

>>> import pywt
>>> x = [3, 7, 1, 1, -2, 5, 4, 6]
>>> db1 = pywt.Wavelet('db1')
>>> cA3, cD3, cD2, cD1 = pywt.wavedec(x, db1)
>>> print(cA3)
[ 8.83883476]
>>> print(cD3)
[-0.35355339]
>>> print(cD2)
[ 4.  -3.5]
>>> print(cD1)
[-2.82842712  0.         -4.94974747 -1.41421356]


If sampling rate is high does that increase the error or decrease the error? Suppose we are tracking a computer mouse on the screen, if we increase the frequency of the sampling rate(to very high values) then would that increase the error in getting the exact position of the mouse or decrease it?

Would we get lots of outliers and we would have to average out?

What I am trying to ask will it get smooth lines or every little flicker in movement will be noted?

Well, i don't see how is that related to DSP or DIP, except of course the issue of sampling, but here is what is going on:

If you are tracking the mouse by latching on an OS' provided "onMouseMove" signal then you are already enjoying (almost) the maximum sampling resolution offered by the combination of your (hardware + OS). That would be the absolute max in Fs accessible by other programs running on the system.

Is that enough? It is, if you are free drawing with small movements trying to get the contour of an object in some photo editing software (for example). But if you click on the mouse and go crazy with it you might start to notice the trick behind the "smooth lines" which is to compose them of linear segments. Therefore, because of the large actual movement compared to the Fs you will start noticing some long lines appearing as the linear segments connect points that are far appart.

Can this Fs be changed? (Increased / Decreased) Yes, it probably can but that would depend on your hardware + software combination. For example, here are a few pieces of information regarding the PS2 mice where you can see a specific command that can be sent to it to change the sampling resolution (search for "set mouse sample rate"). How you actually get to send this command to the mouse so that the sample rate is changed across the system, depends on your OS).

Will higher Fs lead to improved tracking of the actual mouse movement (smaller error)? Yes, that is the point of increasing the Fs, you can track the continuous physical phenomenon much better.

Why don't we crank the Fs up all the way then? Because depending on your Hardware + OS there will be a limit that you can increase the Fs before it starts becoming "noticable" by your system. What you want to avoid is having the OS triggering the "onMouseMove" signal / message so fast that it floods the message queue and at the extremes the system ends up absorbed to serving just this kind of requests.

Similarly, if you are tracking the mouse position with any other means (a camera looking at the screen for "an arrow that points up left" for example, or the same thing by screen grabbing) and possibly with another $Fs_{Track} < Fs_{Sys}$ then the same thing applies. However, in this case, MAYBE you can somehow increase $Fs_{Track}$ by improving the scanning process...but the best case scenario is that you will be approaching $Fs_{Sys}$ not exceeding it.

I hope this helps.

P.S. Regarding "every flicker", modern mice will report a resolution in DPI (Older mice had mechanical rotary encoders and just good enough resolution). Therefore, you can consider your mouse reporting its position over a very fine canvas. In general, mice are good in discriminating movement and will report robust positions over this fine canvas. There are exceptions however where an optical mouse, at rest, will report changes in the position depending on the surface it rests on. In this case, the mouse gets "tricked" and you can avoid this flicker by changing the mat it rests on. For any other means of tracking, the "flicker" will depend on the tracker and Signal to Noise Ratio (SNR).


Real impulse response It will be great if someone can explain me what exactly means "real impulse response".
Further more , what is the effect of such a response on ROC (Laplace plane) and in particular if its restrict the possible area of the poles and zeros?
Thanks

real impulse response

An impulse response that is real-valued.

effect of such a response on ROC

nothing specific, since realness doesn't imply any kind of stability/convergence, far as I can tell. It does tell you something about poles and zeros always appearing in conjugate pairs, but that's it.

restrict the possible area of the poles and zeros

no


Quantization error of SAR-ADCs My question is related to the quantization error and the functionality of SAR-ADCs. In general, the quantization error of an ADC is defined as +-0.5LSB. If this concept is correct a voltage that equals 1.75LSB, would produce 2LSB at the output of the ADC and a voltage that equals 1.25LSB would produce 1LSB at the output of the ADC. Based on the linked example transformation curve we got U_LSB = 1V/(2^3bit) = 0,125V, so an input voltage of 0,21875V (1.75LSB) would produce an output voltage of 0,250V and an input voltage of 0,15625V would produce an output voltage of 0,125V. So far so good. Let’s jump to the basic concept of an SAR-ADC.

An input voltage Vin is sampled by a sample-and-hold block. This sampled voltage is then compared by a comparator with a known voltage. This known voltage V_DAC is provided by the DAC that is supplied by a reference voltage Vref. The digital input for the DAC is generated by a SAR register as a digital word. This digital word depends on the result from the comparator. Initially the SAR begins with setting the MSB. After that the output voltage V_DAC of the ADC is compared with the input voltage Vin. If Vin is greater than V_DAC the MSB stays set. If Vin is less than V_DAC the MSB will be set to zero again. In the next step the SAR register sets the following bit and the comparison process of V_DAC and Vin goes on, until the LSB is reached.
My problem lies in this strict comparison process, which will never allow the output voltage of the ADC to be higher than the input voltage Vin. When I take the example from the first paragraph a SAR-ADC would give me 0,125V (1LSB) as an output voltage, as long as I don’t exceed 0,250V (2LSB). So an input voltage of 1.25LSB and 1.75LSB would both provide 1LSB output. This would mean the concept of +-0.5LSB quantization error wouldn’t be correct for SAR-ADCs and the ideal transformation steps would be shifted to the right for like 0.5ULSB (see diagram with red curve). The new quantization error would be max. -1LSB.

Is this thinking process correct or am I missing something? I can’t imagine that the general concept of quantization error does not match with the architecture of SAR-ADCs…

What the OP shows is the difference between truncation and rounding which in terms of error simply introduces an offset. They generate the same quantization noise spectrally; the standard deviation of the quantization noise (and equivalently noise power) is equivalents since DC offset is ignored in those computations. If there was concern with the offset error introduced (1/2 lsb) it can also easily be subtracted.
Further it should be clarified that DAC's and ADC's will have a finite DC offset regardless due to the analog portion of the electronics, this is often consistent with other DC offsets in the analog chain (consider the offset voltage between the negative and positive Op-amp terminal for example, or carrier leakage in zero-IF receivers), and when an actual issue (for a system where we are trying to measure a DC value) would be incorporated in an overall DC-offset removal algorithm with either calibration or cancellation loop) as part of all offsets combined. Since most DAC's and ADC's have a precision out to what there self-noise factors allow, the ENOB (equivalent number of bits) is often 1 to 2 bits worst than the stated precision (for example, a 12 bit converter may have an ENOB of 10.5 bits). For this reason it is often not even a concern if there is a DC offset introduced of 0.5 bits since that is well below the actual achieved precision overall. These details with regards to ENOB are clarified on ADC and DAC datasheets.
The OP asked for how the noise equivalence can be verified. It is obvious by inspection of the plots given that the only difference between the two is a static offset, but the equivalence can be computed simply enough for further verification by creating a floating point full scale  sine wave and then creating two quantization signals by both rounding and truncating the sine wave to the closest integer. Subtract the quantized fixed-point (integer) signal from the original floating point signal to get the quantization error alone. Then compute the standard deviation of each and they will be equivalent.
This chart I have demonstrates this graphically; specifically what is shown is generalized to a complex waveform (two ADC's with I and Q) but the concept is the same for a single ADC (real waveform), in that the rounded or truncated waveform is the sum of quantization noise with the higher precision waveform (such as floating point or more specifically samples of analog waveform prior quantization through rounding/truncation).

Further, if you use an incommensurate sine wave (choosing a frequency for the sine wave that isn't an integer sub-multiple of the sampling clock), you can also see through subsequent FFT and other statistical analysis how the quantization noise is well modeled as a uniform white noise distribution, and you can also confirm for either case (truncation or rounding) the formula used to estimate the SNR (the total power of the quantization noise floor relative to the full-scale sine wave) as given by:
$$SNR = 6.02 \text{dB/bit} + 1.76 \text{dB}$$
Which I detail further at this link:
What are advantages of having higher sampling rate of a signal?


Addition of signal variance and noise variance Can I add two variances directly? I have my noise variance $\sigma_{noise}=-10 \ \mathrm{dB}$ and signal variance, $\sigma_{signal}=40 \ \mathrm{dB}$. So, can I say the total variance is $\sigma_{total}=30 \ \mathrm{dB}$?

I will use this total variance to generate random variables using a normal distribution with $\mu =0$ and $\sigma = \sigma_{total}$.

If the signal and noise are uncorrelated then the variance of the sum of the two equals the sum of their variances. So, yes, you can add variances (of uncorrelated signals). However, you cannot add dB-values! First of all, a variance in dB does not make sense because dB quantifies the ratio between two values. What you probably mean is dBm or dBW. So before adding variances you have to compute their actual (linear) values from the given values in dB(m/W).


Digital RF signal vector resize algorithm I'm searching for an algorithm for resizing a RF sampled signal in order to draw it without losing any important maximum or minimum values.

The values of the vector are negative/positive around zero.

Probably not what you want, but using the convex hull of your points will preserve the minima and maxima.

R Code Below

#31780
T <- 1000
x <- rnorm(T,0,1)

ix <- seq(1,T)

ix_hpts <- chull(x = ix, y= x)

plot(ix_hpts, x[ix_hpts], type='l', col='red')
points(ix, x)


Sliding window minimum and maximum filter in MATLAB I'm trying to apply a sliding window minimum and maximum filter to an image of a certain window size. Actually, I'm trying to find the optimum window size for it. But I really haven't gotten the hang of it. I presume that I should be using blockproc to implement the sliding window, but not really sure how to find the maximum and minimum filter. As to the implementation itself, should I use loops to slide the window across the entire area of the image ?

ordfilt2 will do this for you:

N=5;
filtered_img = ordfilt2(img, 1, true(N));

for minimum and

filtered_img = ordfilt2(img, N*N, true(N));

for the maximum.  You can also use imdilate and imerode to perform maximum and minimum filters.


What Are the Limitations of a Canny Edge Detector? Broadly a majority of the literature on edge detection algorithms and applications that uses edge detection, references Canny's edge detector.  So much so that it looks like almost "the solution" to edge detection. Certainly, it would do the best job balancing noise and preserving edges.

However, as a simple curiosity, is there an area of concern for Canny's edge detector? or are there areas of applications where Canny will not be best?

In this context, faster implementation is not really concern. The focus of edge detector being good or bad should be the quality and utility of edges generated.

Also, I am really not focusing on implementation specific issues. I am looking for more theoretical limitations or characteristics inherent in the algorithm.

From my experience, the following points are limitations:

The result is binary. You sometimes need a measure of 'how much' the edge qualifies as an edge (e.g. intensity image coming from a Sobel amplitude edge detector)
The amount of parameters leads to infinitely tweaking for getting just that little better result.
You still need to connect the resulting edges to extract the complete edges that seem so obvious for the human eye+mind.
Also due to the gaussian smoothing: the location of the edges might be off, depending on the size of the gaussian kernel.
The method has problems with corners and junctions:

The gaussian smoothing blurs them out, making them harder to detect (same goes for the edges themself)
The corner pixels look in the wrong directions for their neighbors, leaving open ended edges, and missing junctions

This last problem is addressed by the SUSAN method, which connects edges better and also results in nice junctions, as shown by these example figures as given in the linked paper:

Test input image:

Results SUSAN:

Results Canny:

You can clearly see SUSAN finds the corners and junctions instead of Canny.


Circular Time Delay in signal, effect on phase spectra of DFT If you circularly shift a signal x[n]= [-3 -2 -1 0 1 2 3 2 1 0 -1 -2 ] to the right by M=1 compared to M=2, the phase spectra has a lot many zeros when M=2 compared to M=1.

here, I am talking about Circular Time shifting.

I know that if x[n] is even, then DFT is purely real and hence phase is zero. However, why are there way many more zeros when M=2 compared to M=1?

Thanks

If by phase spectrum you mean the sequence $$\bigr(\angle X[0], \angle X[1], \angle X[2], \ldots, \angle X[N-1]\bigr),$$ then remember that a circular shift
by $M$ in the time domain results in each $X[n]$ being multiplied by $\exp(-j2\pi Mn/N)$, and thus causes a change in $\angle X[n]$; in particular, $\angle X[n]$
decreases by $2\pi Mn/N$.  For the case $M=1$, the multipliers are distinct $N$-th roots of unity, that is, the changes in phase are

$$\Bigr(0, 2\pi\frac{1}{N}, 2\pi\frac{1}{N}, \ldots, 2\pi\frac{N-1}{N}\Bigr)$$

all of which are different numbers.  Thus, each entry in the phase
spectrum changes (except for the entry for $n=0$). For the case $M=2$ and $N=12$, the multipliers are $\exp(-j2\pi 2n/12) = \exp(-j2\pi n/6)$ and so the changes in phase are

$$\Bigr(0, 2\pi\frac{1}{6}, 2\pi\frac{2}{6}, \ldots, 2\pi\frac{5}{6},
2\pi\frac{6}{6}, 2\pi\frac{7}{6}, \ldots 2\pi\frac{11}{6}\Bigr)$$

But a phase change of $2\pi\frac{6}{6} = 2\pi$ is the same as
no phase change, a phase change of $2\pi\frac{7}{6}$ is the same
as a phase change of $2\pi\frac{1}{6}$, and so on.  So the displayed
vector above of phase changes is actually

$$\Bigr(0, 2\pi\frac{1}{6}, 2\pi\frac{2}{6}, \ldots, 2\pi\frac{5}{6},
0, 2\pi\frac{1}{6}, \ldots 2\pi\frac{5}{6}\Bigr)$$

which is two periods of a sequence of period $6$.  Note
each entry in the phase spectrum changes except for $n=0$ and
$n=6$ which do  not change.

Now suppose that each $X[n]$ is a real number and thus the phase spectrum
is $0$ everywhere. Then, if you circularly shift by one place ($M=1$ in the
time domain, all entries in the phase spectrum will be nonzero
except for the $n=0$ entry.  If instead you circularly shift by
two places ($M=2$) in the time domain, the phase spectrum will have zeroes in
both the $n=0$ and $n=6$ entries.  Note also that $\exp(-j\pi = -1$,
and so for $M=1$, $X[0]$ and $X[6]$ will still be real-valued (though
there will have been a phase change of $\pi$ in $X[6]$) and similarly, for
$M=2$, $X[0], X[3], X[6]$, and $X[9]$ will still remain real-valued
though there will have been a phase change of $\pi$ in $X[3]$ and $X[9]$.


What is used to initialise difference equations that require data from the output? Given some difference equation for a filter like

$y[k] = ax[k] + bx[k-1]+ cy[k-1]$,

how would you initialize it? Since it needs an old output value (feedback) to calculate the new output value, it seems like the equation will keep on chasing its own tail. I would be tempted to just let

$y[k-1] = x[k-1]$

initially to calculate a value for $y[k]$ after which we could just use that for $y[k-1]$. There might be a better technique for this problem or even a really obvious solution to it, but I am at a loss on how. Most of the digital filter books/papers that I have read have done a great job on explaining how the filter gets derived, but the practical implementation requirements seem to be absent.

This is also my first time posting on Stack Exchange, so any comments on etiquette or style would be appreciated. Just be gentle ;)

If there is no good reason to choose otherwise, you would initialize with zeros, so in your case $y[-1]=0$. Other initializations may be useful in certain situations, e.g., when processing blocks of data to avoid transients between blocks.

Note that if the filter is stable (which it has to be in almost all useful applications), the influence of the initial condition on the output becomes negligible after a while.


How to choose suitable segmentation process for image processing Hello I'm new on image processing.

I know there are a lot of threshold and region growing algorithms to use.

But I would like to know what is the trick for choosing suitable method for an image ?

Currently I'm trying to extract characters from this image

and most of the methods I've tried didn't really satisfy me. So I thought maybe I am doing it wrong. Here are some samples.

Normal Threshold

Sobel

Adaptive Threshold

The thresholding techniques that you would find here should be ideally suited to OCR applications. You could read and easily implement. They are also very efficient.


Is there any reason that cross correlation would perform well or poorly on sparse binary arrays? I am using matlab xcorr to correlate simulated photon count data that has some Gaussian random noise set on top of it and it is working fine when the average value in the arrays is greater than one but when I split the arrays into enough bins such that most elements are zero (conserving the total photon count) and the average is much less than one photon per bin, suddenly the lags from xcorr get much larger (relative to the total length even). For example, at one resolution the lag may be 10 while if I increase resolution ten times, the lags become ~500, much more than ten times the previous value. When I do this without any noise added it works very well, and I am thinking maybe because when the values are mostly zeros and ones and the rare two, any differences such as the xth bin being a one rather than a zero is a large percentage difference in the value and this leads to a large difference in the sums (~integrals) that xcorr calculates. Does anyone have any more knowledge or thoughts on whether this is correct or not?

If you upsample your array to a resolution more than the time uncertainty, or measurement time quantization plus jitter noise, by just zero-filling between events in your event count data, you are adding high-frequency jitter noise that can interfere with correlation results.  You could try upsampling your arrays by using a Gaussian distribution of the appropriate width for each event, instead of just zero-padding between events.  Or, if you have a better statistical distribution for the measurement time offsets/errors, use that.


period of this signal it could be a donkey question but I'm a little confused.

I have this signal $w(t)=sin({\pi}0.1f)$. I have to calculate the period of this signal.

$w(t)=sin(2{\pi}\frac{1}{20}f)$

the period of this signal is: $20s$ while the frequency is $\frac{1}{20}Hz$

is it correct?

You are mostly correct.  The equation should be $w(t) = sin(2\pi\frac{1}{20}t)$.  The variable in the sin should be $t$, not $f$.  The frequency is determined by the $\frac{1}{20}$, not by the variable.

Other than that you are correct.  The period is 20s, which makes it a $\frac{1}{20}$ Hz sin wave.


Sample-rate, filtering, digital-filtering and aliasing I am strugling with a question that I hope someone can help me with.

I am recording single molecule events which I detect is picoampere square deflections.

I wish to use as gentle low-pass bessel filtering as possible.

The lowest filter settings my amplifier allow are 10 kHz and 100 kHz, and my digitizer have a maximal sampling rate of 500 kHz. I am afraid of corrupting my signal to much, but do not have the intuitive understanding of sampling and filtering to know if I am doing something wrong. Here is what I do:

I filter the signal with a 100 kHz bessel filter and digitize it with a 500 kHz sampling rate.
I then wish to filter my digitized data with a 35 kHz digital filter.

Would this mess up my data?
I hear people say that I am on safe ground if i sample at appropximatly 10x my filter settings, but I get to this 'safe zone' only when I do the post-sampling digital filtering. So I guess what I realy do not understand is if the order of filtering, sampling, filtering does something nasty to the data.

I hope I was able to communicate my question clear enough.

Thank you very much,
Best regards,
Michael

An ideal square pulse - which I assume is a model for your up-and-down deflections - would have an infinite bandwidth, but the bulk of the energy is within a bandwidth of $1/T$, where $T$ is the pulse duration.  Roughly speaking, therefore, the 100 kHz Bessel filter will thus allow you to detect pulses of duration $10 \mu s$ and above.  It will also limit the sharpness of the up-and-down transition.

You want the Bessel filter to prevent aliasing, which means that it needs to attenuate signals in the 500 kHz +/- 100 kHz range, which would alias back into the passband of the filter.  By my calculations, the attenuation of a 4th order Bessel filter would be about 23 dB in this range (but check your datasheet), which isn't great but may be good enough in some applications.  It depends on how high of a signal-to-noise ratio you need, how noisy your signal is within this range, and how much spectral content your deflections have in this range (which will depend on their width and the sharpness of their transitions).

The rule-of-thumb you mentioned about sampling 10x above your filter settings may be specific to a particular filter.  In general, you only need to sample at twice the bandwidth of your signal, which is the well-known Nyquist Criteria.  In practice, we usually sample at higher rates because it allows us to use a more realistic, lower cost filter - such as your 4th order Bessel filter.  The required sampling rate, the filter specifications, and the spectral properties of the signal you are sampling are all related.

Having sampled at 500 kHz, you would want to filter your signal with a 35 kHz filter if you are only interested in spectral content less than 35 kHz.  This also allows you to reduce your sample rate to some value greater than 70 kHz, which could reduce the computational burden of whatever processing you are doing.  But you would not be concerned with the 10x rule-of-thumb anymore, because any aliasing due to the original sampling has already occurred.  You can implement a digital filter that is considerably tighter than the Bessel filter.


Multi Resolution STFT for audio I am looking at trying to achieve a multi resolution STFT analysis on a log sine sweep audio measurement.

I have an impulse, and I would like to perform two different windowed length  FFTs on the impulse.

The shorter window length will determine the crossover frequency of the two results. So the this short windowed signal FFT will supply the complex result for x frequency and above.

The longer window will provide all of the results below x.  Although, I am trying to understand how to properly achieve this and synthesise the end result.  Do I need to low pass the longer window before the FFT,  if so, how do I achieve this?

And then how do I put the two results together?

I've done multi-resolution spectrograms before, so they're not impossible (a unicorn).  One method is to zero-pad the shorter data window (for better time locality resolution) to the same length FFT as the longer data window (for greater frequency resolution).  Then you can simply scale for window length, and cross-fade the two sets of FFT magnitude results over some frequency range.

Another possibility (instead of zero-padding) is to use interpolation of the lower frequency FFT or lower time resolution FFT results to create more possibilities for a higher density of spectrogram plot points.

One academic paper suggested computing the local vertical and horizontal contrasts at each spectrogram point and dynamically picking the FFT set with the highest contrast (over the vertical for narrow frequency spectra or across the horizontal for short time events) to plot at each point.

If you keep both (or more!) sets of (overlapped) STFTs, your transform isn't lossy, but redundant, and the resulting image can use that redundancy to obscure less information, usually due to graphical blurring/plotting with typical linear interpolation and/or non-linear color mapping.

Even for reconstruction you can use both sets of data, for instance feeding the short IFFTs to a tweeter and the longer IFFTs to a sub-woofer, with a cross-over in the middle (be careful of phase).


What does "noise power is assumed to be the same at all nodes I was reading a paper about communications between nodes, and in the simulation section they gave the condition of the simulation as follow :

"...The noise power is assumed to be the same at all nodes..."

In the simulation it is never stated that the transmit power of all nodes are not the same, but if you put them so that they are equal you get an error. My guess is that they mean the noise variances $σ^2$ of the nodes are the same. However the simulation result is way off the mark from the result in the paper.

Did I understand their meaning properly?

Although noise is often called channel noise in various block diagrams of radio-frequency (RF) communications systems that communications analysts draw, the truth of the matter is that most actual RF channels are mostly noiseless, and what is referred to as the channel noise is primarily thermal noise in the resistive elements in the front end of the receiver.  The noise power in the receiver that communications engineers have to contend with is thus determined by the temperature of those resistors (and the value of those resistors).  Thus, assuming that the receivers in the various "nodes" are identical and also in the same geographical area so that all the receivers can be expected to be at the same temperature, it is reasonable to use a model in which the noise power is the same in all of the receivers.